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Abstract 

The  major  goal  of  this  research  was  to  investigate  machine  recognition  of  faces.  The 
approach  taken  to  achieve  this  goal  was  to  investigate  the  use  of  the  Karhunen-Loeve  Trans¬ 
form  (KLT)  by  implementing  flexible  and  practical  code.  The  KLT  utilizes  the  eigenvectors 
of  the  covariance  matrix  as  a  basis  set.  Faces  were  projected  onto  the  eigenvectors,  called 
eigenfaces,  and  the  resulting  projection  coefficients  were  used  as  features.  Face  recognition 
accuracies  for  the  KLT  i.oefRcients  were  superior  to  Fourier  based  techmques.  Addition¬ 
ally,  this  thesis  demonstrated  the  image  compression  and  reconstruction  capabilities  of  the 
KLT.  This  thesis  also  developed  the  use  of  the  KLT  as  a  facial  feature  detector.  The  abil¬ 
ity  to  differentiate  between  facial  features  provides  a  computer  communications  interface 
for  non-vocal  people  with  cerebral  palsy.  Lastly,  this  thesis  developed  a  KLT  based  axis 
system  for  laser  scanner  data  of  human  heads.  The  scanner  data  axis  system  provides  the 
anthropometric  community  a  more  precise  method  of  fitting  custom  helmets. 


X 


Face  Recognition  with  the  Karhunen-Loeve  lYansform 


1.  Introduction 


General  Issue 

Face  recognition  is  an  active  area  of  research  at  the  Air  Force  Institute  of  Technology 
(AFIT).  The  AFIT  interest  in  face  recognition  stems  from  the  Department  of  Defense 
(DoD)  need  to  provide  support  for  human  security  efforts.  For  example,  by  scanning  the 
faces  of  people  entering  a  secure  or  restricted  area,  a  face  recognition  system  could  be 
used  to  control  entry.  Similar  defense  related  applications  of  face  recognition  can  easily  be 
found  in  antiterrrorism  and  antinarcotic  operations.  As  an  answer  to  the  DoD  need,  AFIT 
developed  the  Autonomous  Face  Recognition  Machine  (AFRM)  which  scans  a  room  and 
identifies  all  the  faces  in  the  room.  Unfortunately,  the  AFRM  has  limited  accuracy.  In 
order  to  improve  the  accuracy  of  machine  based  face  recognition,  this  thesis  investigates 
the  use  of  the  Karhunen-Loeve  Transform  (KLT). 

The  following  section  provides  the  thesis  background  and  justification.  The  third 
section  formally  states  the  research  problem.  In  support  of  the  problem  statement,  the 
following  sections  are  presented:  research  objectives,  research  questions,  methodology, 
standards,  and  scope.  The  last  section  provides  an  overview  of  the  thesis. 

Background 

The  AFRM  is  the  first  working  face  recognition  system  developed  at  AFIT.  Originally 
completed  in  1985  by  Russel,  subsequent  additions  were  performed  by  Smith,  Lambert, 
Sander,  and  Robb  (18,  21,  10,  19,  14).  To  recognize  subjects,  the  AFRM  windows  the 
face  into  six  regions;  the  AFRM  then  performs  a  Discrete  Fourier  Transform  (DFT)  on 
each  window  (10:3-33-3-36).  The  subject  is  recognized  from  a  set  of  known  images  using 
a  minimum  distance  classification  scheme  (10:3-40).  The  performance  of  the  AFRM,  or 
any  other  pattern  recognition  system,  is  primarily  dependent  on  the  features  used  for 
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classification.  Depending  on  the  test  conditions,  the  AFRM  accuracy  can  vary  from  56% 
to  90%  (14:35-36).  The  poor  accuracy  of  AFRM  can  be  attributed  to  the  use  of  Fourier 
coefficients  as  features.  In  this  instance,  the  Fourier  coefficients  are  not  an  adequate  feature 
set  for  recognition.  It  has  been  shown  that  the  KLT  is  an  optimal  transform  based  on 
minimum  mean  square  error  of  reconstruction  and  maximum  entropy  of  the  representation 
(4:125).  This  thesis  determines  the  optimality  of  the  KLT  feature  set  for  face  recognition. 

Problem  Statement 

The  primary  problem  is  to  recognize  human  faces.  The  key  task  associated  with  this 
problem  is  determining  a  good  set  of  features  that  would  allow  a  machine  to  accurately 
recognize  human  faces. 

Research  Objectives 

The  objective  of  this  research  is  to  investigate  the  optimality  of  the  KLT.  The  sec¬ 
ondary  research  objective  is  to  develop  an  architecture  and  methodology  for  face  recogni¬ 
tion  with  the  KLT. 

Research  Questions 

The  research  completed  for  this  thesis  answers  the  foUowing  questions: 

1.  Can  the  KLT  be  implemented  in  a  numerically  practical  manner? 

2.  Does  the  Karhunen-Loeve  orthogonal  transformation  of  the  face  provide  a  better 
feature  set  for  face  recognition  than  traditional  Fourier  techniques? 

3.  What  recognition  accuracies  can  be  expected  from  a  Karhunen-Loeve  basis  set? 

4.  Does  holistic  face  recognition,  which  utilizes  the  whole  face  as  a  feature,  provide 
better  results  than  localized  AFRM  techniques? 

5.  Can  a  KLT  based  system  be  used  as  a  computer  communications  interface  for  the 
non-vocal? 

6.  Can  the  KLT  be  used  as  an  axis  system  for  the  anthropometry  community? 
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Methodology 

As  paxt  of  this  thesis,  face  recognition  software  is  developed  on  the  Silicon  Graphics 
4D  Personal  Iris  workstation.  The  system  software,  which  is  written  in  ANSI  C,  is  devel¬ 
oped  using  a  modular  design  approach.  The  algorithms  are  sufficiently  flexible  to  be  used 
on  faces  as  well  as  any  other  image  recognition  problem,  e.g.  tanks,  trucks,  and  jeeps. 

This  research  investigates  the  KLT  and  develops  software  to  implement  the  KLT 
and  Fourier  Transform.  This  thesis  also  demonstrates  the  image  coding  and  reconstruc¬ 
tion  capabiflties  of  the  KLT.  Most  importantly,  the  KLT  coefficients  will  be  used  for  face 
recognition  and  the  result  will  be  compared  to  more  traditional  Fourier  techniques.  The 
recognition  will  utilize  an  Euclidean  distance  or  first  nearest  neighbor  classifier.  Lastly, 
two  additional  applications  of  the  KLT  will  be  investigated.  The  first  KLT  application 
investigates  the  abilities  of  the  KLT  to  differentiate  facial  features  which  would  permit  the 
KLT  to  be  used  as  a  facial  feature  communications  interface  for  the  non-vocal.  The  second 
application  investigates  the  correlation  properties  of  the  Karhunen-Loeve  for  use  as  an  axis 
system  for  the  anthropometry  community. 

Standards 

The  most  important  performance  criteria  is  classification  accuracy  which  is  the  per¬ 
centage  of  correct  classifications  out  of  the  total  number  of  inputs.  The  inputs  consist  of 
both  known  and  unknown  subjects. 

Scope 

The  focus  of  the  thesis  is  feature  selection  and  recognition.  This  thesis  assumes  that 
the  detection  and  segmentation  phases  of  the  recognition  problem  are  adequately  met  by 
the  AFRM.  The  assumption  is  valid  if  the  whole  face  is  used  as  a  feature  and  Lambert’s 
moving  target  indicator  is  used. 


1-3 


Overview  of  Thesis 


This  chapter  provides  a  brief  background  and  defines  the  problem,  objectives,  ques¬ 
tions,  methodology,  standards,  and  scope  of  the  research.  Chapter  2  provides  additional 
background  and  a  review  of  current  research.  The  methodology  used  to  solve  the  problem 
is  presented  in  Chapter  3.  The  results  of  the  research  are  presented  in  Chapter  4.  Chapter 
5  summarizes  the  significant  results  and  provides  recommendations.  Appendix  A  lists  the 
source  code  utilized  for  the  thesis. 
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//.  Background 


Introduction 

This  chapter  provides  background  and  evaluates  past  and  current  research  in  the  area 
of  face  recognition.  Early  face  recognition  systems  approached  the  face  recognition  problem 
from  an  individual  facial  feature  perspective,  for  example,  the  Fast  Fourier  Transform 
(FFT)  of  the  eyes  or  the  distance  between  the  eyes  (14,  5).  Those  early  systems  suffered 
from  poor  accuracy  and  segmentation.  Recent  research  has  approached  face  recognition 
from  a  more  holistic  approach  by  using  the  whole  face  as  a  feature  (27,  2). 

This  chapter  begins  with  an  evaluation  of  the  Air  Force  Institute  of  Technology 
(AFIT)  Autonomous  Face  Recognition  Machine  (AFRM)  (14).  To  emphasize  that  a  good 
face  recognizer  requires  good  features,  L.  D.  Harmon’s  profile  recognizer  is  evaluated  (5). 
Lastly,  the  chapter  concludes  with  the  evaluation  of  two  more  recent  holistic  face  recogni¬ 
tion  techniques.  At  the  Massachusetts  Institute  of  Technology  (MIT),  Turk  and  Pentland 
have  developed  an  approach  using  a  holistic  Karhunen-Loeve  Transform  (KLT)  (27).  The 
second  approach  is  from  the  University  of  California  at  San  Diego,  where  Cottrell  has 
developed  a  neural  network  based  face  recognition  technique  (2). 

Air  Force  Institute  of  Technology  Autonomous  Face  Recognition  Machine 

The  design  goal  of  the  AFIT  AFRM  is  to  reliably  recognize  faces  without  human 
intervention.  Face  recognition  which  is  a  relatively  simple  task  for  humans  can  not  be 
reliably  accomplished  by  machines  without  a  great  deal  of  human  intervention.  The  AFRM, 
AFIT’s  first  attempt  at  face  recognition,  began  with  Routh’s  Ph’d  dissertation  on  the 
Cortical  Thought  Theory  (CTT).  The  CTT  states  that  the  brain  uses  the  center  of  mass 
calculation,  or  gestalt,  of  an  im^e  for  identification  and  storage  (16).  With  Routh’s 
theoretical  work  in  place,  the  AFRM  progressed  from  initial  construction,  by  Russel,  to 
subsequent  modifications  by  various  AFIT  students  (14, 10). 

The  AFRM,  which  consists  of  a  video  camera  and  a  computer,  first  acquires  the 
image.  Acquiring  the  ims^e  consists  of  placing  the  subject  in  the  field  of  view  of  the 
camera.  Segmentation,  or  locating  the  face  in  the  image,  is  accomplished  with  one  of 
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two  algorithms.  The  Smith  algorithm  searches  throughout  the  image  for  the  brightness 
‘signature’  of  a  face  (14:4).  For  instance,  the  brightness  signature  of  a  face  is  the  change 
in  brightness  that  occurs  as  the  camera  pans  across  the  eyes.  The  second  acquisition 
algorithm  uses  a  moving  target  indicator,  also  called  a  spatio-temporal  filter,  to  segment 
the  moving  pixels  of  the  face  from  the  stationary  background  (10:3-2).  With  acquisition 
and  segmentation  complete,  preprocessing  must  be  accomplished  to  provide  normalized 
images  for  feature  extraction  and  recognition. 

AFRM  preprocessing  consists  of  brightness  normalization,  contrast  enhancement, 
and  scaling.  One  of  the  most  significant  results  of  Lambert’s  work  is  the  development 
of  a  brightness  normalization  technique  referred  to  as  ‘Lambertization’  (10:E-2  -  E-3). 
‘Lambertization’  consists  of  scanning  an  image  pixel  by  pixel  and  computing  the  local 
average  of  a  pixel  and  its  neighbors.  The  localized  average  is  subtracted  from  that  pixel 
and  its  neighbors.  The  resulting  image  is  normalized  with  respect  to  local  brightness.  With 
the  image  normalized,  the  AFRM  begins  feature  extraction  and  recognition. 

To  extract  features,  a  window  is  drawn  around  the  face.  The  window  is  divided  into 
six  sub-regions.  A  gestalt  is  calculated  on  each  of  the  six  sub-regions  (10:2-25).  A  later 
modification,  successfully  made  by  Robb,  substitutes  the  low  frequency  Fourier  coefficients 
for  the  gestalt.  This  substitution  is  valid  if  you  consider  Mahler’s  animal  cracker  experi¬ 
ment,  which  validated  Kabrisky’s  theory  that  the  low  frequency  Fourier  coefficients  are  a 
good  model  for  the  human  visual  recognition  system  (12,  6).  The  gestalt  or  Fourier  coef¬ 
ficients  are  then  stored  in  the  data  base  for  all  known  subjects.  With  features  extracted, 
the  remaining  task  is  recognition  and  classification. 

To  recognize,  the  AFRM  acquires  and  segments  a  test  face  until  the  gestalt  or  Fourier 
coefficients  are  determined.  The  Euclidean  distance  between  the  test  subject  and  every 
subject  in  the  data  base  is  calculated.  A  minimum  distance  implies  a  match  between  the 
test  and  data  base  subjects. 

The  performance  of  the  AFRM  can  be  best  measured  by  Robb’s  results  in  1989.  Using 
uncontrolled  lighting  and  background,  Robb’s  accuracy  varied  between  56%  and  90%  with 
a  test  set  of  50  subjects  (14:35-36).  Although  the  AFRM  performance  is  satisfactory, 


2-2 


the  accuracy  is  insufficient  for  many  real  applications.  The  following  techniques  seek  to 
dramatically  improve  recognition  performance  by  finding  a  better  set  of  features. 

Profiles  as  Features 

Can  a  set  of  features  be  found  that  permits  accurate  recognition  of  faces?  In  L.D. 
Harmon’s  Machine  Identification  of  Faces,  this  question  is  answered  (5).  Harmon  utilizes 
facial  profiles  as  features  for  face  recognition.  Figure  2.1  exhibits  the  angles  and  distances 
between  landmarks  used  by  Harmon  as  features.  Using  principal-component  analysis,  also 
called  eigen  or  covariance  analysis,  Harmon  isolates  15  optimal  features  out  of  38  features 
(5:108).  Harmon’s  recognition  accuracy  is  on  the  order  of  98-100%,  with  a  population  of 
112  manually  segmented  subjects  (5:104).  Harmon’s  results  prove  that  with  good  features, 
an  accurate  face  recognizer  can  be  built.  Harmon’s  work  is  being  reevaluated  by  Hiroshi 
Aga,wa  as  a  means  of  coding  three  dimensional  face  images  (1). 

Holistic  Face  Recognition 

The  AFRM  and  Harmon’s  profile  recognizer  are  traditional  feature  based  systems. 
Each  part  of  the  face  is  a  feature,  and  as  a  result,  a  larger  amount  of  measurement  error  is 
introduced.  In  holistic  face  recognition,  the  whole  face  is  a  feature  which  in  turn  results  in 
less  measurement  error  and  improved  recognition  accuracy.  For  example,  instead  of  using 
interpupillary  distance  as  a  feature,  holistic  features  would  consist  of  the  Fourier  coefficients 
of  the  entire  face.  There  are  two  current  techniques  that  implement  holistic  approaches  to 
face  recognition.  At  MIT,  Pentland  and  Turk  generate  orthogonal  eigenfaces  with  the  KLT 
(27).  At  UCSD,  Cottrell  generates  nonorthogonal  basis  functions  with  a  neural  network 
(2). 

Karhunen-Lo^ve  Transforms 

Eigenvalues.  Eigenvalue  analysis  is  a  very  powerful  engineering  tool.  In  image  pro¬ 
cessing,  eigen  analysis  supports  image  coding  and  compression  applications.  Eigenvalues 


2-3 


Figure  2.1.  Hannon’s  profile  with  critical  locations  for  feature  extractions  annotated. 
(5:97) 
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are  the  solution  for  A  in  the  equation 


CxX  =  xX  (2.1) 

where  Cx  is  a  by  matrix,  £  is  an  eigenvector,  and  A  is  a  scalar  called  the  eigenvalue 
(9:345).  The  set  of  non-trivial  solutions  to  the  above  equation  are  the  eigenvalues  which 
are  also  referred  to  as  the  spectrum  of  Cx  (9:346).  Each  eigenvalue  has  a  corresponding 
eigenvector,  x.  The  eigenvalues  are  found  from  the  solutions  to  the  following  determinant 

|Cx-AI|£=0  (2.2) 

where  I  is  the  identity  matrix  (9:347).  The  resulting  eigenvalues  are  used  to  calculate  the 
eigenvectors.  As  a  by-product  of  the  eigenvalue  determination,  the  resulting  eigenvectors 
are  orthogonal. 

Theory  of  the  Karhunen-Lo^ve  Transform.  The  KLT  is  an  orthogonal  trans¬ 
form  similar  to  the  Fourier  Transform.  Unlike  Fourier,  which  has  a  basis  set  consisting  of 
sines  and  cosines,  the  KLT  basis  set  consists  of  the  eigenvectors  of  the  covariance  matrix. 
The  KLT  is  also  guaranteed  to  provide  the  most  efficient  and  accurate  transformation  by 
minimizing  mean  square  error  (MSE)  in  reconstruction  and  maximizing  entropy  of  the 
representation. 

The  KLT  can  be  calculated  in  the  three  basic  steps.  First,  obtain  the  mean  and 
covariance  function  from  the  images.  Second,  calculate  the  eigenvalues  for  the  centered 
covariance  matrix,  where  centered  implies  that  the  mean  has  been  subtracted  off  the  image. 
For  each  eigenvalue,  calculate  the  respective  eigenvector  and  rearrange  the  eigenvectors 
into  descending  eigenvalue  order.  Select  only  the  k  largest  eigenvalued  eigenvectors  out  of 
all  eigenvectors.  This  has  the  effect  of  greatly  reducing  the  dimensionality  of  the  pattern 
recognition  problem,  i.e.  simplifying  an  n*  problem  to  k,  while  minimizing  MSE.  The 
resulting  k  eigenvectors  are  the  new  KLT  basis  set.  Lastly,  the  image  can  be  transformed 
into  the  new  basis  set  (3). 
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Before  be^nning  the  KLT  explanation,  a  definition  is  in  order.  First,  define  the  n 
by  n  pixel  image  as  a  vector 

where  x'  is  the  vector  describing  the  image  and  x\  is  the  first  pixel  of  the  image. 
The  construction  of  x*  is  shown  in  Figure  2.2.  The  mean  vector,  fn^,  is  calculated  by 


Figure  2.2.  Definition  of  image  vector,  x',  where  the  first  pixel  is  x^  and  the  last  pixel 


averaging  each  pixel  over  all  M  faces  in  the  training  set 

M 

TTlx  =  1/M^X*  (2.4) 

where  M  is  the  number  of  faces  in  the  training  set  and  x*  is  the  vector  describing  the 
face  (4:123).  The  construction  of  m*  is  shown  in  Figure  2.3. 
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where  x\  is  the  first  pixel  of  the  image.  If  the  mean  image  is  subtracted  off  all  of  the 
images,  then 


/  _t  _t  _t  _« 


Cx  =  1/M 


Em  »  I 
»=1 


(2.7) 


The  covariance  matrix  represents  the  information  content  or  relative  change  in  each  image. 
By  finding  the  covariance  between  the  first  pixel  and  every  other  pixel  in  the  image,  a 
judgement  can  be  made  as  to  the  amount  of  information  in  the  first  pixel.  For  example, 
a  pixel  that  does  not  change  relative  to  other  pixels  provides  little  information,  while 
a  pixel  with  high  covariance  provides  more  information.  Finding  the  eigenvectors  and 
eigenvalues  of  the  covariance  matrix  transforms  the  covariance  matrix  into  an  orthogonal 
space.  In  effect,  the  relative  pixel  changes  are  mapped  into  an  orthogonal  space.  To  an 
engineer,  this  orthogonal  transformation  represents  an  optimum  distance  between  vectors 
and  minimum  decision  errors  which  is  the  ultimate  pattern  recognition  goal.  The  last  KLT 
step  of  selecting  the  k  highest  valued  eigenvectors  creates  a  basis  set  with  only  the  largest 
variance  eigenvectors. 

The  eigenvalue  solution  to  the  covariance  matrix  is  an  simultaneous  equation. 
This  fairly  lengthy  computation  can  be  accomplished  with  an  approximation  which  wiU  be 
demonstrated  in  a  future  section. 

As  previously  mentioned,  the  resulting  eigenvalues  and  eigenvectors  are  rank  or¬ 
dered  with  the  highest  eigenvalued  vector  in  ro'v  one.  Generally,  the  number  of  eigenvectors 
to  be  kept  for  reconstruction  or  recognition  is  determined  using  the  MSE.  The  MSE  be¬ 
tween  the  reconstruction  x  and  the  original  im^e,  z,  is 

n*  k 

=  (2.8) 

i=i  >=i 

which  simplifies  to 

n* 

MSE=  Y,  (2-9) 
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where  Xj  is  the  eigenvalue,  n?  is  the  total  number  of  eigenvalues,  and  k  is  the  number 
of  eigenvectors  used  for  reconstruction  (4:125).  Using  the  MSE,  the  k  largest  eigenvalued 
eigenvectors  are  selected  resulting  in 


^  eii  ci2  ...  ein»  ^ 


u  = 


*21  C22 


*211* 


\  etl  Cfc2  •••  / 


(2.10) 


where  u  is  the  truncated  KLT  matrix,  e/h,  is  the  value  of  the  eigenvector,  and 
the  rows  represent  the  eigenvectors  associated  with  the  Xk  eigenvalue  (4:124).  Each  row 
represents  a  single  eigenface,  and  k  basis  functions  define  the  space.  In  vector  form,  the 
original  image  x  is  transformed  by 


y/t  =  ttfc(f  -  mx)  k=l...M  (“2. 11) 

where  y*  is  the  element  of  the  projection  vector,  y,  in  KLT  space,  is  the  fc*** 
eigenvector,  x  is  the  original  image,  and  mg  is  the  mean  face  (4:125).  In  KLT  space, 
y  represents  an  image  in  a  coordinate  system  defined  by  the  eigenvectors.  To  reconstruct 
the  image,  multiply  y  by  the  inverse  transformation,  u~^  =  u^,  and  add  the  mean  face 

m 

i  =  (^  Ufc^y*)  +  mx  (2.12) 

k=l 

X  is  the  reconstructed  image,  Uk  is  the  k*^  eigenface,  yk  is  the  k*^  KLT  coefficient,  rux  is 
the  average  face,  and  M  is  the  number  of  images  in  the  trsuning  set.  With  a  theoretical 
understanding  of  the  KLT,  an  application  by  Turk  and  Pentland  can  be  evaluated. 

Turk/Pentland  Eigenface.  The  Matthew  A.  Turk  and  Alex  P.  Pentland  eigenface 
is  an  attempt  to  recognize  faces  as  a  holistic  entity  (27:43).  The  Turk/Pentland  approach 
is  heuristically  simple.  First,  with  the  KLT,  develop  an  optimal  basis  set  for  recognition 
(27:45).  With  the  basis  set,  transform  the  images  and  develop  a  data  base  with  the  known 
subjects  and  their  KLT  coefficients. 
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To  recognize  an  unknown  subject,  tbe  test  image  is  transformed  into  KLT  space, 
and  the  Euclidean  distance  of  the  KLT  coefficients  to  each  known  subject  is  calculated.  If 
the  minimum  distance  meets  a  certain  threshold,  the  image  is  classified  as  known.  If  the 
image  does  not  meet  the  minimum  distance  criteria,  the  face  is  classified  as  unknown.  For 
example.  Figure  2.4  shows  a  two  dimensional,  two  eigenvectors  (ui,U2),  face  space, 
which  falls  into  class  Qi  is  classified  as  known.  On  the  other  hand,  yj  which  does  not  fall 
into  the  circled  re^on  of  a  known  class  is  classified  as  unknown.  Lastly,  represents  a 
non-face  (26:49)  because  it  is  relatively  far  from  face  space.  This  type  of  classification  can 
be  readily  implemented  with  a  nearest  neighbor  classifier  or  with  a  more  robust  Adaptive 
Resonance  Theory  (ART)  Neural  Network  (15).  With  that  explanation,  a  more  analytical 
evaluation  of  Turk  and  Pentland  follows. 


_ ^2 _ 

Figure  2.4.  A  simplified  2-D  face  space  where  Ui  and  U2  are  the  eigenface  basis  set  and 
the  known  individuals  are  (27:49) 


KLT  Approximation.  The  KLT  is  numerically  cumbersome.  According  to  Turk, 
the  by  solution  of  the  covariance  matrix  can  be  greatly  simplified  if  the  number  of 
points  in  the  image  space,  i.e.  n^,  is  much  greater  than  the  M  training  faces  (26:6).  For 
example,  a  training  set  of  40  subjects,  each  128  by  128  pixels,  would  satisfy  the  condition 
that  M  «  where  M  =  40  and  =  128^.  The  following  technique  is  based  on  the 
previous  condition. 
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Turk  and  Pentland  simplify  the  KLT  by  estimating  Cx-  First,  if  A  is  defined  as  the 
matrix  which  contains  each  training  face  as  a  column 

A=($i|  $2|  ...  (2.13) 

where  is  the  statistically  centered  training  face.  The  first  statistically  centered 
training  face  would  equal  the  first  face  minus  the  mean  face,  .  Typically,  the 

KLT  is  created  by  first  calculating 


Cx  =  AA***  (2.14) 

where  the  eigenvectors  of  Cx  are  determined  using  Equations  2.1  and  2.2.  Substituting 
AA^  for  Cx  in  Equation  2.2  yields 


AA'^^u*  =  UkX  (2.15) 

where  A  is  the  training  matrix,  A^  is  the  transpose  matrix,  ti*  is  the  eigenvector,  and  A 
is  the  eigenvalue.  In  effect,  Uk  represents  a  vector  that  when  multiplied  by  the  covariance 
matrix,  Cx  =  AA"^,  equals  the  same  eigenvector  scaled  by  A. 

Now  consider  the  eigenvectors  of  L  =  A'^A,  the  resulting  matrix  is  only  M  by  M. 
An  eigen  solution  to  A'^A  would  result  in 

A'^Aui  =  vifii  (2.16) 


where  the  eigenvectors  of  L  are  and  the  eigenvalues  are  pj.  If  both  sides  aire  premultiplied 

by  A,  the  result,  which  is  very  reminiscent  of  the  solution  to  Equation  2.15,  is 


AA^AiJ; 


Avifii 


(2.17) 


where  the  eigenvectors  of  AA*^  are  Avj-  and  the  eigenvalues  are  m.  Now  substitute  Cx 
for  i4i4^ 


CxAiTj  =  AviHi 


(2.18) 
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Ai^'  and  m  are  the  respective  eigenvectors  and  dgenvalues  of  the  covariance  matrix  Cx- 
Therefore,  finding  vj,  the  eigenvectors  of  L,  and  multiplying  vi  by  A  results  in  Uk,  the 
desired  eigenvectors  of  Cx  (26:7). 

To  summarize,  first  generate  A  with  the  training  faces.  Second,  multiply  A*^  and 
A  to  get  L.  Next,  find  the  eigem'alues  and  eigenvectors  of  L,  and  multiply  the  resulting 
eigenvectors,  t^,  by  A  to  get  the  desired  eigenvectors,  Uk,  of  the  original  covariance  matrix 

M 

lifc  =  (2.19) 

i=i 

where  is  the  training  face,  Vkj  is  the  element  of  the  eigenvector  of  L,  and  M  is 
the  number  of  training  faces.  For  instance,  the  first  eigenface,  tti,  is  a  linear  combination  of 
each  element  in  the  first  eigenvector,  lii,  with  its  corresponding  training  face  in  A  (26:6-7). 
The  matrix  L  has  simplified  the  eigen  analysis  from  to  M  (27:7). 

With  the  u  transformation  matrix,  Turk  transforms  the  image  into  face  space  with 
the  following  transformation 


Un  =  Unix*  -  mx)  n  =  1  •  •  •  (2.20) 

where  x*  is  the  image,  rux  is  the  mean  face,  u^  is  the  eigenface,  u;„  is  the  coefficient 
of  the  n*^  eigenface,  and  k  is  the  total  number  of  truncated  elgenfaces  (27:46).  The  values 
of  Un  determine  the  location  of  the  transformed  image  in  multidimensional  face  space. 
Lastly,  Turk  calculates  the  Euclidean  distance  to  determine  which  face  class  is  nearest  the 
unknown  face 

€j  =  |f2  -  Qjl  (2.21) 

where  Cj  is  the  distance  to  each  known  face  class  ,  Q  is  the  unknown  face  class,  and  Qj  is 
the  vector  for  the  known  face  class  (27:47).  Pentland  first  classifies  an  unknown  face 
as  either  face  or  not  face,  based  upon  cy  being  below  a  specified  threshold.  The  minimum 
(j  determines  class  membership. 

The  Pentland  and  Turk  system  utilizes  a  Sun  workstation.  The  working  system 
generates  a  seven  dimensional  eigenface  face  space.  The  recognition  accuracy  is  96%  for  a 
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population  of  16  subjects  (27:19).  Figure  2.5  displays  Pentland’s  reconstruction. 


Figure  2.5.  On  the  left  is  the  original  image.  On  the  right  the  Turk  and  Pcntland  recon¬ 
structed  image.  (27:47) 


According  to  Pentland,  the  approach  does  have  limitations.  Pentland's  KLT  is  sensi¬ 
tive  to  head  size,  which  forces  size  normalization.  A  limitation  in  segmentation  also  e.xisls 
requiring  Pentland  and  Turk  to  use  a  spatio-temporal  filter  similar  to  the  moving  target 
indicator  of  the  AFRM  (27:52-53). 

Neural  Networks  and  Face  Recognition 

Garrison  W.  Cottrell,  at  UCSD,  implements  a  neural  network  to  recognize  face.s 
(2).  Cottrell’s  neural  network  is  another  holistic  face  recognition  technique.  Cottrell 
first  trains  his  neural  network  using  a  back  propagation  neural  network  configured  as  an 
autoassociator,  Figure  2.6  a.  The  input  into  the  autoassociator  is  a  61  by  61  pixel  image. 
At  the  output,  the  autoassociator  monitors  the  MSE  between  the  input  and  the  outiiut. 
The  autoassociator  iteratively  adjusts  the  weights  of  the  hidden  layer  nodes  until  the  M.SE 
is  minimized.  Minimizing  M.SE  results  in  the  input  being  reproduced  at  the  output.  Once 
the  MSFi  is  minimized,  training  is  complete.  The  outputs  of  the  hidden  nodes  are  fi'aturc's 
that  are  used  for  recognition  by  another  neural  network. 
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To  recognize  faces,  the  input  layer  and  hidden  layer  of  the  autoassociator,  Figure 
2.6  b,  are  input  into  a  backpropagation  network  (BPN).  The  BPN  is  trained  to  recognize 
the  nonorthogonal  projections  of  each  subject.  If  a  test  face  which  is  input  into  the  two 
layer  network,  Figure  2.6  b,  activates  one  of  the  outputs  of  the  BPN,  the  test  subject 
is  classified  as  known.  If  the  test  face  does  not  activate  a  BPN  output,  the  test  image  is 
classified  as  unknown. 

The  limitations  of  Cottrell’s  network  are  similar  to  Pentland’s.  However,  Cottrell’s 
neural  network  approach  has  one  key  advantage  over  the  eigenface  approach.  The  neural 
network  technique  does  not  have  the  computational  burden  of  the  KLT  (27:53).  This 
advantage  is  eliminated  if  the  KLT  is  implemented  on  a  neural  network  as  demonstrated 
by  Tarr  (24). 


Figure  2.6.  After  the  above  Cottrell  autoassociator  is  trmned,  the  autoassociator  input 
and  hidden  layers  are  used  as  the  input  to  a  second  network  which  in  turn 
classifies  the  faces. 
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Conclusion 


This  chapter  evaluated  AFIT’s  AFRM,  Harmon’s  face  profile  research,  Pentland  and 
Turk’s  eigenfaces,  and  Cottrell’s  autoassociator.  All  the  research  implies  that  the  pri¬ 
mary  factor  in  recognition  accuracy  is  feature  selection.  More  recent  research  indicates 
that  holistic  face  recognition  provides  improved  results  over  individual  facial  feature  tech¬ 
niques.  With  these  considerations,  this  thesis  wiU  concentrate  on  exploiting  holistic  KLT 
approaches  to  face  recognition. 
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III.  Methodology 


General  Issue 

One  of  the  objectives  of  this  thesis  is  to  evaluate  the  Karhunen-Loeve  Transform 
(KLT)  as  a  possible  basis  set  for  face  recognition.  The  evaluation  of  the  KLT  is  divided 
into  three  phases.  The  first  phase  consists  of  using  the  KLT  with  only  one  training  image 
presented  multiple  times.  The  purpose  of  the  first  phase  is  to  simplify  the  problem  in 
order  to  evaluate  the  algorithmic  implementation  of  the  KLT  in  C.  The  second  phase 
consists  of  image  reconstruction  using  both  small  and  large  population  training  sets.  The 
third  phase  utilizes  the  KLT  for  recognition  and  classification.  As  a  by-product  of  the 
three  phases,  an  evaluation  of  finding  or  segmenting  faces  in  KLT  space  will  be  performed. 
Additionally,  a  facial  feature  communicator  for  the  non- vocal  and  a  KLT  based  axis  system 
will  be  implemented.  Two  additional  KLT  applications  will  be  evaluated.  Before  the  three 
phases  are  presented,  this  chapter  be^ns  with  an  overview  of  the  primary  software  routines 
utilized  for  the  thesis  effort. 

Code  development 

Introduction.  The  software  written  in  ANSI  C  consists  of  three  main  algorithms. 
The  first  algorithm,  kl-transform2.c,  calculates  the  KLT  of  the  input  training  set.  With  the 
eigenfaces  generated  by  kl-transform2.c,  the  second  algorithm,  recon.c,  calculates  the  KLT 
coefficients  and  reconstructs  the  input  image.  The  last  algorithm,  findjnin.c,  classifies 
the  coefficients  based  on  Euclidean  distance.  In  addition  to  the  previous  routines,  two 
KLT  preprocessing  routines,  gwind.c  and  CJateS.c,  will  be  detailed.  The  software  routine 
gwind.c  windows  the  images  with  a  Gaussian.  The  second  preprocessing  routine,  C  Jate5.c, 
centers  all  the  images  relative  to  a  reference.  AH  software  listings  for  this  thesis  can  be 
found  in  Appendix  A. 

KLT  software.  The  primary  function  of  the  KLT  program,  kLtransform2.c,  is  to 
calculate  the  optimum  basis  set,  or  eigenvectors,  and  the  average  vector  of  a  training  set. 
The  trsuning  set  used  to  create  the  eigenvectors  defines  the  basis  set,  therefore  a  training 
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set  of  faces  generates  an  optimum  basis  of  eigenfaces.  Note  that  the  eigenfaces  are  only 
optimum  for  images  in  the  training  set.  As  the  images  become  less  like  the  training  images, 
the  optimality  quickly  disappears.  This  implies  that  a  basis  set  of  eigenfaces  would  only 
be  suitable  for  faces.  If  a  KLT  basis  set  is  desired  for  tanks,  for  example,  the  training  set 
would  have  to  be  modified  to  include  tanks.  The  resulting  basis  set  would  then  provide 
eigentanks.  The  routine  kl.transform2.c  is  primarily  used  on  faces  in  this  thesis,  but  the 
software  can  readily  be  implemented  with  any  training  set  of  images. 

The  KLT  software  is  developed  with  two  key  references:  Turk  and  Pentland’s  Eigen¬ 
faces  for  Face  Recognition  asid  Gonzalez  and  Wintz  Digital  Image  Processing  (26,  4).  The 
C  program  kl-transform2.c  which  is  documented  in  Appendix  A  is  executed  from  the  Unix 
command  line  with 

kl.transform  tralningfile  sizeofimage  numberoftrain 

where  tralningfile  is  the  name  of  the  file  which  lists  the  faces  in  the  KLT  training  set, 
sizeofimage  is  the  total  image  size  in  pixels,  and  numberoftrain  is  the  number  of  images  in 
the  training  set.  The  block  diagram  of  kl-transform2.c  is  shown  in  Figure  3.1. 

First,  each  training  face  has  the  average  face  subtracted  off  resulting  in 

A=(r,|  Fal  ...  IFm)  (3.1) 

where  A  is  the  KLT  training  set  matrix,  Fm  is  the  statistically  centered  image,  and 
M  represents  the  total  number  of  training  images.  To  simplify  the  calculation  of  the 
KLT,  rather  than  calculating  the  eigenvectors  of  the  covariance  matrix  Cx  =  AA^,  the 
eigenvectors  of  L  are  found 

L  =  A'^A  (3.2) 

The  Numerical  Recipes  in  C  routine  jacobi.c  finds  the  solution  to  Equation  2.2,  where 
L  is  substituted  for  Cx  (13).  A  second  Numerical  Recipes  in  C  routine  eigsort.c  sorts 
the  eigenvectors  generated  by  jacobi.c  in  decreasing  eigenvalue  order  (13).  Note  that  the 
resulting  eigenvectors  from  L  are  only  of  length  M  which  is  the  size  of  the  training  set. 
One  additional  step  is  required  to  arrive  at  the  appropriate  length  eigenvectors.  This 
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is  accomplished  with  Equation  2.19  to  find  uic  the  length  eigenvector.  The  outputs  of 
kl_transfonn2.c  are  the  average  face,  eigenfaces,  and  eigenvalues. 


training  facet 


avg_face.dat 

elgenfac801.dat 


Figure  3.1.  Block  diagram  for  the  C  program  kl-transform2.c  which  finds  the  eigenfaces, 
i.e.  KLT  basis  set,  and  the  averse  face  for  a  given  training  set. 


Reconstruction  of  faces.  Thus  far  the  eigenfaces  and  average  face  have  been 
calculated  with  the  program  kl-transform2.c.  This  subsection  explains  the  implementation 
of  the  software  routine  recon.c.  The  program  recon.c  calculates  the  KLT  coefRcients  of 
the  image,  and  with  the  eigenspace  representation  of  the  image,  reconstructs  the  original 
image. 

The  C  program  recon.c  which  is  documented  in  Appendix  A  is  invoked  at  the  Unix 
prompt  as  follows 

rscon  input inags  numberinrecon  sizaofinags 
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where  inputimage.gra  is  the  image  to  be  reconstructed,  inputimage.rec  is  the  reconstruc¬ 
tion,  numberinrecon  is  the  number  of  eigenfaces  used  in  the  reconstruction,  and  sizeofimage 
is  the  total  number  of  pixels  in  the  input  image.  First,  recon.c  reads  the  input  image,  the 
average  face,  and  the  eigenfaces  needed  for  reconstruction.  The  image  x  is  then  projected 
into  the  new  eigen  coordinate  system.  The  new  vector  y  is  created  by 

yfc  =  ui(x  -  m*)  k  =  1...M  (3.3) 

where  yt  is  the  element  or  coefficient  of  the  k*^  eigenface,  tii  is  the  k*'*  eigenface,  x  is 
the  input  image,  mj.  is  the  average  face,  and  M  is  the  number  of  images  in  the  training 
set.  The  vector  y  is  the  projection  of  the  input  image  into  KLT  space.  These  coefficients 
are  also  features  for  recognition  and  classification.  The  reconstruction  is  the  weighted  sum 
of  the  coefficients  of  y  with  the  appropriate  eigenface.  In  vector  form 

m 

i  =  (H  ^Jvk)  +  rfix  (3.4) 

k=t 

where  x  is  the  reconstructed  image,  Uk  is  the  k*^  eigenface,  yk  is  the  k*'*  KLT  coefficient, 
ntx  is  the  average  face,  and  M  is  the  number  of  images  in  the  training  set.  The  outputs  of 
recon.c  are  the  reconstructed  image  and  the  KLT  coefficients. 

Euclidean  Distance  Classifier.  The  previous  two  software  routines  provided  an 
eigenface  basis  set,  generated  KLT  coefficients,  and  reconstructed  an  image  from  its  KLT 
coefficients.  To  actually  recognize  test  from  prototype  faces,  a  classifier  must  be  imple¬ 
mented.  In  this  thesis  a  Euclidean  distance  or  first  nearest  neighbor  (1-NN)  classifier  is 
utilized.  The  routine  find_min.c  documented  in  Appendix  A  is  executed  from  the  Unix 
command  line  with 

find.ain  trainfila  tastfila  t.train  i.tast  f.faaturas 

where  trainfile  is  a  file  containing  the  prototype  vectors,  testfile  is  a  file  containing  the  test 
vectors,  #.train  is  the  total  number  of  prototype  vectors,  #.test  is  the  total  number  of  test 
vectors,  and  #  Jeatures  is  the  number  of  features  in  each  vector.  Each  test  and  prototype 
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bnags 

1 

Caleulato 

Wsights 

^  •igonfaes.dat 

avgLfac«.dat 

1 

Rsconstruct 

KLT  coaffleisnts  to  classifier 

1 

rseonstnictsd  Imags 

Figure  3.2.  Block  diagram  for  C  program  recon.c 


vector  ends  with  a  string  or  label  identifying  the  vector  class.  The  program  findjnin.c 
statistically  normalizes  the  test  and  prototype  features.  The  normalized  Euclidean  dis¬ 
tance  from  the  test  vector  to  each  prototype  vector  is  calculated.  The  minimum  distance 
from  test  to  prototype  implies  that  the  test  vector  is  of  the  same  class  as  the  prototype. 
After  program  completion,  find_min.c  lists  each  test  vector  with  the  associated  winning 
prototype  and  distance.  The  probability  of  error  of  a  Euclidean  distance  or  1-NN  classifier 
approximates  a  Bayesian  classifier  as  the  number  of  test  and  prototype  vectors  approaches 
infinity  (25:83).  In  practice,  with  more  than  twelve  classes  the  1-NN  approximates  the 
Bayesian  classifier. 

Preprocessing  Paces  with  a  Gaussian  Window.  Applying  a  Gaussian  window 
to  face  images  has  several  advantages.  A  Gaussian  windowed  face.  Figure  3.3,  has  its  center 
accentuated  and  the  outline  of  the  face  de-emphasized.  From  a  recognition  prospective, 
de-emphasizing  the  outline  of  the  head  de-emphasizes  hair  and  background  and  emphasizes 
the  center  of  the  face.  Both  effects  are  desireable  for  recognition.  From  a  similar  biological 
perspective,  when  humans  recognize  faces,  the  eyes  spend  more  time  at  the  center  of  the 
face  and  less  on  the  outline(ll).  Figure  3.4  demonstrates  a  typical  eye  scan  pattern  of  a 


3-5 


human  examining  a  face.  The  eye  scan  pattern  which  is  taken  with  an  occulometer  is  in  a 
sense  Gaussian  with  considerably  more  scans  on  the  center  of  the  face  than  on  the  outline. 
An  additional  benefit  of  Gaussian  windowing  is  that  the  centering  of  the  images  is  focused 
on  the  eyes  instead  of  the  outline  of  the  head.  In  short,  Gaussian  windowing  the  face  is 
biologically  motivated,  provides  more  desirable  recognition,  and  more  desirable  centering. 


Figure  3.3.  An  example  of  a  face  windowed  with  a  Gaussian. 

The  window  is  implemented  with  the  software  routine  gwind.c,  see  Appendix  A. 
The  program  takes  an  input  image  and  multipUes  each  pixel  by  a  corresponding  pixel  in 
the  Gaussian  window,  Figure  3.5.  First  the  routine  gwind.c  prompts  for  the  image  size, 
and  the  routine  then  prompts  for  the  x  mean,  x  variance,  y  mean,  and  y  variance.  The 
center  of  the  Gaussian  window  is  defined  by  x  mean  and  y  mean  while  x  variance  and 
y  variance  define  the  amount  of  change  in  the  x  and  y  direction.  For  example.  Figure 
3.5  demonstrates  a  Gaussian  window  with  parameters:  x  mean=65,  x  variance=35,  y 
mean=55,  y  variance=50.  The  center  of  the  128  by  128  pixel  Gaussian  window  is  right  65 
pixels  and  up  55  pixels,  and  the  window  variances  define  the  tightness  of  the  window.  In 
most  instances,  the  mean  values  stay  constant  and  the  variances  are  varied. 

Preprocessing  Faces  by  Centering.  The  KLT  is  not  shift  invariant.  In  other 
words,  shifting  a  face  in  a  scene  provides  different  KLT  coefficients.  This  undesirable 
effect  can  be  eliminated  by  centering  the  images  relative  to  a  reference.  The  centering 
operation  is  performed  by  a  simple  routine  CJateS,  see  Appendix  A.  The  program  does  a 
Fast  Fourier  Transform  correlation  of  the  input  image  and  an  already  centered  reference 
face.  The  correlation  indicates  how  much  shift  is  necessary  in  x  and  y  pixels  to  align  the 
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Figure  3.4.  The  ori^nal  scene  of  the  little  ^rl  is  on  the  left  and  the  resulting  occulometer 
pattern  is  on  the  right.  Note  the  concentration  of  scans  on  the  center  of  the 
face.  (11) 


input  image  with  the  reference.  The  routine  then  shifts  the  input  image  the  appropriate 
amount  and  stores  the  shifted  image.  Unless  otherwise  indicated,  the  program  CJateS 
does  not  implement  energy  normaliztion. 

Testing.  To  test  the  code,  a  simple  KLT  problem  is  devised.  The  training  set  con¬ 
sists  of  nine  instantiations  of  a  single  face,  with  each  face  differing  by  the  addition  of  various 
amounts  of  Gaussian  noise.  Properly  running  software  provides  accurate  reconstruction 
and  a  rapidly  decaying  mean  square  error  (MSE)  plot.  The  results  of  this  test  are  provided 
in  the  next  chapter. 
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(a)  (b) 


Figure  3.5.  a.  A  Gaussian  window  with  parameters  x  mean=65,  x  variance=35,  y 
mean=55,  and  y  variance=50.  b.  A  Gaussian  window  with  parameters  x 
mean=65,  x  variance=25,  y  mean=55,  and  y  variance=30. 

KLT  Reconstruction 

Small  Population  Reconstruction  using  the  KLT.  The  KLT  is  typically  seen 
in  text  books  as  an  image  compression  technique  (4).  To  evaluate  the  KLT,  a  test  is  set  up 
with  a  small  population  of  six  faces  in  the  training  set.  The  faces  are  recorded  on  video  tape 
and  digitized  using  Snapshot,  a  NeXT  application  package.  The  input  images,  320  by  480 
pixels,  are  scaled  to  256  by  256  pixels.  To  prevent  blurred  reconstructions,  all  the  images 
are  centered.  A  Gaussian  window  is  then  applied  which  decreases  the  effects  of  hair  line, 
head  size,  and  background  on  recognition.  A  second  centering  operation  is  performed  to 
concentrate  centering  on  the  center  of  the  face  ,  i.e.  eyes.  With  all  the  images  aligned,  the 
KLT  is  implemented  with  kl-transform2.c.  As  previously  outlined,  kl.transform2.c  creates 
the  average  face  and  eigenfaces.  After  the  eigenfaces  and  average  face  are  calculated,  any 
of  the  six  faces  can  be  reconstructed  with  recon.c.  The  results  of  the  reconstruction  are 
provided  in  the  next  chapter. 

Large  Population  Reconstruction  using  the  KLT.  The  reconstruction  of  im¬ 
ages  from  large  KLT  populations  seems  to  be  a  logical  extension  to  the  small  population 
test.  The  primary  objective  of  the  following  tests  is  to  evaluate  the  quality  of  reconstruction 
with  different  size  training  sets  and  Gaussian  windows. 

The  first  test  evaluates  a  training  set  of  27  subjects.  Each  subject  has  three  different 
instantiations  for  a  total  of  81  faces  in  the  training  set.  The  Gaussian  window  of  Figure 
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3.5  a  is  applied.  This  test  evaluates  reconstruction  of  images  in  and  out  of  the  training  set. 

The  second  test  evaluates  the  same  training  set  as  in  test  one  but  with  a  tighter 
Gaussian  window  around  the  faces.  The  Gaussian  window  of  Figure  3.5  b  is  applied.  This 
test  evaluates  the  effects  of  different  window  sizes  on  reconstruction  quality. 

The  third  test  evaluates  27  trsdning  faces  with  only  one  instantiation  per  test  sub¬ 
ject.  The  Gaussian  window  of  Figure  3.5  a  is  applied.  This  test  evaluates  the  effects  on 
reconstruction  of  using  only  one  version  of  each  training  face. 

The  fourth  test  evaluates  reconstruction  with  only  18  faces  in  the  training  set  and  the 
Gaussian  window  of  Figure  3.5  b  is  applied.  This  test  evaluates  trmning  set  size  reduction 
and  image  reconstruction  quality. 

The  fifth  test  demonstrates  the  effects  of  not  using  Gaussian  windowed  images. 
Lastly,  the  sixth  test  evaluates  reconstruction  of  a  subject  not  in  the  KLT  training  set. 

The  basic  set-up,  for  all  the  tests  is:  develop  the  KLT,  reconstruct  faces,  and,  lastly, 
heuristically  evaluate  the  quality  of  the  reconstruction  for  each  test.  The  results  of  the 
large  population  reconstruction  are  provided  in  Chapter  4. 

KLT  and  Recognitioa 

Introduction.  Previous  experiments  evaluated  the  reconstruction  capabibties  of 
the  KLT,  but  the  goal  of  this  research  is  to  develop  a  set  of  features  for  accurate  face 
recognition. 

Single  Person  Verification.  The  first  evaluation  of  the  KLT  recognition  capability 
consists  of  single  person  verification.  For  example,  a  subject  enters  a  controlled  access 
facUity  and  claims  to  be  person  A.  The  subject  then  enters  into  a  terminal  a  personal 
identification  number  (pin)  code.  At  the  same  time  a  video  camera  takes  a  picture  of 
the  face.  The  KLT  coefficients  of  the  image  verifies  if  person  A  is  truly  person  A,  or  an 
impostor  with  the  right  pin. 

This  problem  is  evaluated  as  a  two  class  problem  with  class  one  consisting  of  the  test 
subject  A  and  class  two  consisting  of  every  other  person  in  the  data  base.  A  KLT  with 
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18  different  training  images  encodes  the  images.  Using  only  the  first  five  KLT  coefficients, 
a  Backpropagation  Network  (BPN)  with  momentum  classifies  the  data  (22).  Various  ex¬ 
periments  are  run  varying  the  number  of  input  features  and  hidden  layer  nodes.  Lastly, 
using  the  same  test  scenario  and  data,  a  Euclidean  distance  classifier  measures  recognition 
accuracy.  The  results  of  single  person  verification  are  provided  in  Chapter  4. 

Multi-person  Recognition.  Determining  the  identity  of  a  subject  from  multiple 
subjects  is  the  multi-person  recognition  problem.  The  population  of  faces  consists  of  fifty- 
five  subjects.  Each  subject  is  diptized  into  128  by  128  pixel  frame  and  is  part  of  the 
KLT  trsuning  set.  All  the  images  are  shifted,  aligned,  and  windowed  as  in  Figure  4.5. 
Executing  kl-transform2.c  provides  the  average  face,  eigenfaces,  and  eigenvalues.  With  the 
eigenfaces,  the  features  can  be  extracted  using  recon.c. 

Before  the  features  are  extracted,  the  six  images  are  divided  into  test  and  prototype 
sets.  The  test  set  consists  of  two  images  which  will  be  compared  to  the  remaining  four 
in  the  prototype  set.  The  prograun  recon.c  extracts  the  features  providing  the  KLT  coef¬ 
ficients  of  each  image.  With  the  KLT  coefficients  of  both  the  test  and  prototype  sets,  a 
Euclidean  distance  classifier,  findjtnin.c,  compares  the  test  set  to  the  prototype  set.  A  min¬ 
imum  distance  implies  a  match.  The  aforementioned  multi-person  test  is  performed  again, 
except  only  forty  images  are  included  in  the  KLT  training  set.  The  results  of  multi-person 
recognition  are  provided  in  Chapter  4. 

As  a  means  of  comparing  the  KLT  to  a  more  common  technique,  the  low  frequency 
FFT  magnitude  coefficients  are  used  for  recognition.  The  FFT  is  energy  normalized  and 
the  zero  frequency  response  (DC)  is  the  center  of  the  FFT.  Since  the  FFT  magnitude  of 
the  face.  Figure  3.6,  is  symmetric,  only  the  top  half  of  the  FFT  magnitude  is  utilized. 
For  example,  a  first  order  block  of  the  FFT  magnitude  consists  of  a  block  with  boundaries 
defined  by  DC  plus  one  coefficient  to  the  right,  left,  and  up.  The  first  order  block  provides 
six  coefficients  because  the  block  is  three  by  two. 


3-10 


Figure  3.6.  FFT  magnitude  of  a  face  with  the  first  order  3  by  2  block 


Finding  Faces  in  Face  Space 

Finding  a  face  in  a  scene  is  the  first  step  in  any  pattern  recognition  system.  The  Air 
Force  Institute  of  Technology  (AFIT)  Autonomous  Face  Recognition  Machine  (AFRM), 
implements  a  spatio  temporal  filter  or  moving  target  indicator  to  segment  a  face  from  the 
overall  scene.  Segmention  based  on  the  KLT  would  segment  on  ‘faceness’  which  would  be 
more  desireable  than  segmentation  based  on  motion.  Two  approachs  will  be  evaluated. 

The  first  segmentation  approach  will  scan  a  scene.  Each  sub-image  scanned  will  be 
projected  into  face  space.  If  the  Euclidean  distance  of  the  KLT  coefficients  of  the  sub¬ 
image  is  less  than  a  certain  threshold,  a  face  is  probably  present.  This  approach  is  the 
same  approach  taken  by  Pentland  (27). 

The  second  segmentation  approach  will  take  the  FFT  of  the  scene  and  of  the  eigen- 
faces.  An  energy  normalized  correlation  between  the  scene  and  eigenface  is  performed  to 
determine  the  location  of  the  face  in  the  scene.  The  results  of  both  techniques  are  provided 
in  Chapter  4. 
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Karhunen-Lo^ve  based  Facial  Feature  Communicator  for  the  Non-vocal 

Many  times  a  Masters  thesis  only  serves  to  enhance  academic  understanding.  In 
this  instance,  however,  an  application  of  the  Karhunen-Lo^ve  may  provide  a  child  with 
cerebral  palsy  a  means  of  communicating  with  the  outside  world.  The  child  has  cerebral 
palsy  and  is  incapable  of  speaking.  Although  the  child  can  not  speak,  she  can  manage 
facial  expressions.  The  following  test  evaluates  the  ability  of  the  KLT  eigenvectors  or 
eigen  images  to  differentiate  between  two  expressions.  The  first  expression,  tongue  out,  is 
classified  as  a  ‘yes’.  The  second  expression,  mouth  open,  is  classified  as  a  ‘no’. 


Two  128  by  128  pixel  images  of  the  subject  shown  in  Figure  3.7  are  selected  as 
prototype  images.  The  two  prototype  images  are  used  as  training  images  to  generate  the 


Figure  3.7.  The  two  prototype  images  of  the  non-vocal  subject.  On  the  left  is  the  tongue 
out  or  ‘yes’.  On  the  right  the  mouth  open  or  ‘no’. 


KLT  eigenvector  basis  set.  The  prototype  images  are  first  Gaussian  windowed  with  Figure 
3.8  to  remove  the  effects  of  the  background  and  enhance  the  mouth  area.  The  images 
are  then  aligned  with  CJateS.c  a  FFT  based  correlation  algorithm.  Finally,  the  KLT  is 
calculated  with  kl.transform2.c  of  Appendix  A  generating  the  average  vector  and  the  first 
two  eigenvectors.  A  block  diagram  of  the  aforementioned  process  is  shown  in  Figure  3.9. 

With  the  KLT  basis,  the  next  step  consists  of  calculating  the  KLT  coefficients.  Figure 
3.10  summarizes  the  following  procedure.  An  input  image  is  first  Gaussian  windowed  with 
Figure  3.10.  The  image  is  then  aligned.  The  KLT  coefficients  are  calculated  with  the 
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Figure  3.8.  The  Gaussian  window  used  on  the  images  of  the  non-vocal  subject. 


Prototype  Images 

* 

Gaussian  Window 

H 

Align  Images 

KLT 

i 

Eigenface 
Average  face 

Figure  3.9.  The  process  used  to  generate  the  KLT  basis  set  from  the  prototype  images 
of  the  non- vocal  subject. 
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routine  from  Appendix  A  recon.c.  Calculating  the  KLT  coefficients  is  done  with  Equation 
3.3  and  consists  of  subtracting  off  the  mean  from  the  input  image  and  then  taking  the  dot 
product  between  the  image  and  each  of  the  two  eigenvectors.  The  two  KLT  coefficients 
are  classified  with  the  Appendix  A  routine  find_min.c,  a  Euclidean  distance  classifier. 
Before  the  test  images  of  Figure  3.11  are  classified,  the  prototype  images  are  processed 
through  Figure  3.10  to  determine  the  prototype  KLT  coefficients.  With  the  prototype 
KLT  coefficients  for  the  ‘yes’  and  ‘no’  facial  expressions,  the  test  images  can  be  projected 
into  Karhunen-Lobve  space  and  classified  with  the  NN  classifier.  The  minimum  distance 
to  a  prototype  determines  if  the  input  or  test  ims^e  is  a  ‘yes’  or  ‘no’.  The  results  of  this 
test  are  presented  in  Chapter  4. 


Test  Images 

Eigenfaces 
Average  face 
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Align  Images  |  — ► 

KLT  feature  extractor  | 

NN  Classifier  | 

Yes  or  No 

Figure  3.10.  The  process  used  to  calculate  the  KLT  coefficients  and  classify  the  images 
as  ‘yes’  or  ‘no’. 


Ksirhunen-Loeve  Axis  System 

The  Karhunen-Loeve  Transform  (KLT)  is  typically  utilized  in  image  compression  or 
recognition  applications.  However,  if  the  KLT  is  to  be  truly  useful  as  a  block  transform, 
it  is  essential  to  be  able  to  correlate  using  the  KLT.  A  sponsor  of  this  research  had  a 
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Figure  3.11.  The  two  top  test  images  are  mouth  open  or  ‘no’.  The  two  bottom  test 
images  are  tongue  out  or  ‘yes’. 
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unique  data  set  of  human  heads  taken  with  a  laser  radar  scanner.  The  scanner  distance 
measurements  of  the  head  are  used  to  fit  flight  gear  on  pilots.  The  laser  scanner  provides 
detailed  measurements  for  all  of  the  features  of  the  face  and  head. 

A  current  short  coming  of  the  laser  scanner  data  collection  is  the  lack  of  an  axis 
system.  An  axis  system  speeds  the  search  of  features  by  placing  all  the  heads  in  the 
same  relative  position.  Like  a  map  with  a  lon^tude  and  latitude  coordinate  system,  the 
axis  system  simplifies  the  task  of  finding  difierent  land  marks  on  the  head.  Therefore, 
searching  for  noses  will  always  occur  in  the  same  general  area  of  the  axis  system.  The 
following  research  demonstrates  the  use  of  Karhunen-Lohve  eigenvectors  or  eigen  images 
as  an  axis  system  and  compares  the  results  to  am  equivalent  Fourier  baised  system. 

Laser  Scan  Data  The  laser  scan  data  of  the  head  consists  of  distance  measurements 
of  the  head.  The  scanner  goes  around  the  head  and  collects  the  data.  The  training  images 
will  define  the  basis  set;  therefore,  selection  of  the  training  set  should  be  representative  of 
the  laiser  scan  data.  The  data  consists  of  a  matrix  which  is  256  data  points  in  the  vertical 
direction  and  512  data  points  in  the  horizontal  direction.  The  distance  measurements  are 
converted  from  distance  to  gray  scale  using  the  routine  kljned.c  which  is  documented  in 
Appendix  A.  The  converted  gray  scale  512  by  256  pixel  image  is  shown  in  Figure  3.12. 

First  a  KLT  training  set  of  nine  laser  scanner  head  images  is  chosen  at  random  to 
represent  the  entire  population.  The  trmning  images  are  aligned  manually  to  improve 
the  quality  of  the  resulting  eigen  images.  The  KLT  is  calculated  using  the  C  program 
kl.transform2.c  which  is  listed  in  Appendix  A.  The  routine  kl-transform2.c  provides  the 
Karhunen-Lohve  eigenvectors  and  average  vector. 

With  the  basis  set  calculated,  a  non-trsuning  image  Figure  3.12  a,  is  selected  and 
shifted.  After  every  shift,  the  projection  into  Karhunen-Loeve  space  is  calculated  using 
only  the  first  eigenvector.  The  projection  is  calculated  with  the  Appendix  A  routine  recon  .c 
which  uses  E^quation  3.3  where  tii  is  the  first  eigen  image,  z*  is  the  shifted  image,  and 
u)\  is  the  KLT  coefficient.  The  goal  of  the  test  is  to  demonstrate  that  a  correlation  peak 
would  exist  when  the  test  image  is  aligned  with  the  eigenvector  or  eigen  image.  The 
head  scan  image  is  shifted  at  20  pixel  increments,  300  pixels  to  the  right  and  to  the  left. 
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Figure  3.12.  Two  laser  scans  of  the  head.  The  top  is  the  original  512  by  256  pixel  image, 
and  the  bottom  is  the  shifted  to  axis  system  image. 
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Simultaneously,  the  image  is  shifted  at  10  pixel  increments,  90  pixels  up  and  down.  Once 
the  peak  is  found,  the  test  image  is  a^ain  shifted  at  one  pixel  increments  to  find  a  precise 
peak. 

As  a  means  of  comparison,  the  FFT  correlation  between  the  head  scanner  image  and 
the  first  eigen  image  is  calculated.  A  correlation  peak  implies  that  the  head  scan  is  aligned 
with  the  eigen  image.  The  results  of  this  test  are  presented  in  Chapter  4. 

Conclusion 

This  chapter  first  presented  various  KLT  and  support  algorithms  developed  for  this 
thesis.  The  chapter  then  presented  the  methodology  used  to  investigate  the  reconstruction 
of  images  using  the  KLT  with  both  large  and  small  population  training  sets.  Thirdly,  the 
chapter  provided  the  methodology  used  for  face  recognition  using  the  KLT  and  FFT.  The 
last  three  sections  presented  the  methodology  used  to  investigate  a  face  finder,  a  facial 
feature  communicator  for  the  non-vocal,  and  a  KLT  based  axis  system.  The  results  of  the 
methodology  are  presented  in  the  next  chapter. 
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IV.  Results 


Code  Development  and  Testing 

In  order  to  test  the  two  main  algorithms,  kl-tran8form2.c  and  recon.c,  the  following 
test  is  devised.  The  first  algorithm,  kl-transform2.c,  performs  the  Karhunen-Loeve  Trans¬ 
form  (KLT).  The  kl-transform2.c  training  set  includes  nine  instantiations  of  the  same  face, 
Figure  4.1  a.  The  nine  faces  differ  only  by  the  addition  of  various  amounts  of  Gaussian 
noise.  After  program  execution,  the  resulting  eigenfaces.  Figure  4.1  b  -  d,  demonstrate 
that  most  of  the  information  or  energy  is  in  the  first  eigenface.  The  mean  square  error 


(b)  (c)  (d) 


Figure  4.1.  (a)  The  KLT  training  set  of  nine  instantiations  of  same  face  with  various 

amounts  of  Gaussian  noise  (b)  first  eigenface  of  the  nine  images  (c)  second 
eigenface  (d)  third  eigenface 

(MSE)  plot  of  Figure  4.2  also  indicates  that  a  major  part  of  the  energy  or  information  is 
in  the  first  eigenvector.  Since  all  of  the  images  are  basically  the  same  face  perturbed  by 
Gaussian  noise,  it  is  reasonable  to  conclude  that  only  one  eigenface  is  needed  to  represent 
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all  ^line  faces.  This  example  confirms  the  operation  of  the  code  and  more  importantly  pro¬ 
vides  insight  into  the  KLT.  The  first  eigenface  represents  most  of  the  face  energy  whereas 
the  second  and  third  eigenfaces  represent  the  orthogonal  projections  of  the  changes  in¬ 
duced  by  the  addition  of  Gaussian  noise  (7).  Figure  4.4  demonstrates  a  two  dimensional 
eigenspace  where  it  can  be  seen  that  most  of  the  energy  in  the  faces  can  be  projected  onto 
the  first  eigenface  and  the  changes  introduced  by  the  Gaussian  noise  can  be  taken  into 
account  by  the  second  orthogonal  eigenface. 


Number  of  Eigenvectors 


Figure  4.2.  Mean  squared  error  of  reconstruction  of  Figure  4.3  a  -  c 


Figure  4.3.  (a)  Reconstruction  using  one  eigenface  (b)  using  two  eigenfaces  (c)  using  all 

nine  eigenfaces 

The  second  algorithm  test  of  recon.c  utilizes  the  eigenfaces  from  kLtransform2.c.  A 
face  is  selected  from  the  nine  training  faces  and  is  provided  as  an  input  to  recon.c.  The 
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Figure  4.4.  An  example  of  a  two  dimensional  eigenspace  where  aU  nine  images  are  the 
same  face  with  various  amounts  of  Gaussian  noise.  The  first  eigenface,  shown 
as  a  bold  vector,  represents  most  of  the  energy  in  the  faces,  and  the  second 
eigenface  accounts  for  the  Gaussian  perturbations.(7) 

program  recon.c  generates  the  KLT  coefficients  and  reconstructs  the  original  image  as 
shown  in  Figure  4.3.  As  indicated  earlier,  the  MSE  plot  of  Figure  4.2  indicates  that  the 
reconstruction  should  only  need  one  eigenface.  Figure  4.3  confirms  that  the  reconstruction 
is  acceptable  using  only  one  eigenface. 

KLT  Reconstruction 

Small  Population  Reconstruction.  The  test  setup.  Figure  4.5,  consists  of  six 
training  faces,  each  scaled  from  a  320  by  480  image  to  256  by  256  pixels.  The  program 
CJate5.c  of  Appendix  A  aligns  the  images.  A  Gaussian  window  is  applied  with  gwind.c 
of  Appendix  A  and  a  second  alignment  with  CJate5.c  concentrates  on  the  eyes,  nose, 
and  mouth  of  the  images.  As  discussed  in  Chapter  2,  the  second  alignment  improves  the 
reconstruction  and  recognition  performance.  Using  kLtransform2.c,  the  KLT  creates  the 
average  face,  eigenvalues,  and  six  eigenfaces. 

Figure  4.6  indicates  10%  MSE  reconstruction  results  when  using  three  eigenfaces. 
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Figure  4.5.  The  KLT  process  consists  of  first  defining  a  training  set,  centering  the  images 
with  CJate5.c,  windowing  with  gwind.c,  centering  once  again,  and  lastly 
performing  the  KLT  with  kl_transform2.c.  The  result  is  the  average  face, 
eigenface,  and  eigenvalues  (not  shown).  Note  all  source  code  is  listed  in 
Appendix  A. 
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(a)  (b)  (c)  (d) 


(e)  (f)  (g) 


Figure  4.7.  Small  population  reconstruction  using  (a)  original  image;  reconstruction  us¬ 
ing  (b)  one  eigenface  (c)  two  eigenfaces  (d)  three  eigenfaces  (e)  four  eigenfaces 
(f)  five  eigenfaces  (g)  six  eigenface 
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An  image  is  first  reconstructed  using  one  eigenface,  Figure  4.7  b,  and  the  result  has  no 
resemblance  to  the  original  of  Figure  4.7  a.  By  the  addition  of  the  third  eigenface,  the 
reconstruction  of  Figure  4.7  d  appears  acceptable.  The  addition  of  eigenfaces  four  through 
six,  Figure  4.7  e-g,  adds  very  little  to  the  quality  of  the  reconstruction. 

The  KLT  permits  the  reconstruction  of  any  of  the  six  faces  using  only  three  of  the 
eigenfaces.  Representing  an  image  by  only  three  coefficients  represents  a  dramatic  com¬ 
pression  of  approximately  21,000:1.  In  a  communications  sense,  instead  of  transmitting  256 
by  256  pixels,  only  three  coefficients  need  to  be  transmitted  provided  that  both  sender  and 
receiver  have  copies  of  the  first  three  eigenfaces  and  average  face.  However,  this  significant 
compression  is  costly.  First,  the  compression  is  limited  to  a  small  population  of  images, 
six  in  th'.s  instance.  Second,  a  significant  amount  of  processing  is  necessary  to  determine 
the  eigenfaces.  Additionally,  both  sender  and  receiver  must  have  copies  of  the  eigenfaces. 
The  pre  and  post  processing  impact  would  be  less  worrisome  if  the  population  could  be 
increased.  Phase  two  of  KLT  reconstruction  consists  of  larger  population  reconstruction 
experiments. 

Large  Population  Reconstruction.  The  first  test  evaluates  a  KLT  training  set  of 
three  128  by  128  pixel  images  per  subject  and  a  population  of  27  subjects.  The  images  are 
Gaussian  windowed  with  Figure  3.5  a  which  has  the  following  parameters:  x  mean=65,  x 
variance=35,  y  mean=55,  y  variance=50.  This  test  evaluates  the  size  of  the  training  set  and 
reconstruction  quality  of  images  in  and  out  of  the  training  set.  The  reconstruction,  using 
20  eigenfaces,  is  shown  in  Figure  4.8.  The  first  column  is  the  original  image.  The  second 
column  is  the  reconstruction  of  a  KLT  training  set  image.  The  third  column  is  an  image  of 
the  subject  not  in  the  KLT  training  set.  The  results  indicate  that  the  KLT  reconstruction 
for  a  large  population  is  much  worse  than  the  small  population  reconstruction.  Column 
three  indicates  that  the  KLT  can  reconstruct  images  not  in  the  training  set  provided  the 
subject  is  in  the  KLT  training  set. 

The  second  test  evaluates  the  effects  of  the  Gaussian  window  on  the  reconstruction. 
The  window  is  tightened  as  in  Figure  3.5  b  to:  x  mean=65,  x  variance=25,  y  mean=55, 
y  variance=30.  The  resulting  quality  is  unaffected  by  window  changes.  Figure  4.9. 
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Figure  4.9.  Test  2:  reconstruction  with  20  eigenfaces  of  large  population,  81  training 
images,  and  tight  Gaussian  window 
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The  third  test  to  improve  reconstruction  quality  varies  the  amount  of  training  faces 
used  per  subject  from  three  to  one,  for  a  total  of  twenty-seven  pictures.  The  faces  are 
Gaussian  windowed  as  in  the  first  test  with:  x  mean=65,  x  variance=35,  y  mean=55, 
y  variance=50.  The  results  of  Figure  4.10  are  similar  to  those  of  the  first  test,  but  a 
close  examination  of  the  images  reveals  a  mar^nai  decrease  in  reconstruction  quality.  For 
example,  note  the  loss  of  detail  around  the  eyes. 

The  fourth  test.  Figure  4.11,  demonstrates  the  use  of  a  smaller  set  of  eighteen  training 
subjects  with  the  following  Gaussian  window:  x  mean=65,  x  vairiance=25,  y  mean=55,  y 
variance=30.  Once  again,  mar^nal  decrease  in  reconstruction  quality  results. 


Original  Training  Non-Training 

Image  Image  Image 


Figure  4.10.  Test  3:  Reconstruction  with  20  cigenfaces  of  large  population,  27  training 
images,  and  only  one  image  per  subject 

The  fifth  test.  Figure  4.12,  demonstrates  the  effects  of  not  using  Gaussian  windowed 
images.  Without  a  Gaussian  window,  the  alignment  of  CJateb.c  performs  poorly.  A  poor 
alignment  results  in  the  blurred  reconstruction  seen  in  Figure  4.12.  This  test  demonstrates 
the  importance  of  centering  or  aligning  the  images,  and  it  proves  the  shift  sensitivity  of 
the  KLT. 
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Figure  4.11.  Test  4:  Reconstruction  with  20  eigenfaces  of  large  population,  18  training 
images,  and  a  tight  Gaussian  window 


Figure  4.12.  Test  5:  Reconstruction  of  non-Gaussian  windowed  images 
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The  sixth  test,  Figure  4.13,  demonstrates  a  reconstruction  of  a  subject  not  in  the 
KLT  training  set.  A  training  set  of  81  images  was  used  to  create  the  KLT  basis  set.  Note 
that  the  reconstruction  has  no  resemblance  to  the  original  image.  This  result  is  somewhat 
disturbing  because  it  indicates  that  the  KLT  did  not  generalize  for  reconstruction.  In  other 
words,  subjects  not  in  the  KLT  training  set  can  not  be  reconstructed. 


Original  Reconstruction 

Image  Non-Training  Subject 


Figure  4.13.  Test  6:  Reconstruction  using  20  dgenfaces  of  subject  not  in  the  KLT  training 
set. 


Overall,  the  best  reconstruction  appears  to  be  that  of  test  one,  a  large  training  set 
with  three  images  per  subject  in  the  KLT  training  set.  In  reality,  the  quality  of  all  four 
reconstructions  is  good  and  is  extremely  better  than  previous  KLT  image  reconstruction  of 
faces  (27).  For  the  first  test,  which  is  the  best  reconstruction,  the  compression  is  16384:20 
or  approximately  800:1  for  81  training  images.  This  is  a  significant  image  compression 
figure.  The  image  compression  comes  at  the  cost  of  constrained  population,  reconstruction 
quality,  and  computational  overhead. 

The  reconstruction  results  are  enlightening  and  valuable  from  an  image  processing 
point  of  view,  but  the  true  goal  of  the  research  is  face  recognition.  A  transform  that  re¬ 
constructs  images  very  well  does  not  necessarily  make  a  good  image  recognition  transform, 
and  vice  versa.  Therefore,  the  most  important  results  of  the  research,  face  recognition  and 
classification,  are  to  follow. 
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KLT  and  Recognition 

Single  Person  Verification.  Single  person  verification  consists  of  recognizing  a 
persons  identity.  For  instance,  if  a  subject  cl^ms  to  be  John  Doe,  does  his  face  match 
John  Doe’s?  This  problem  is  a  two  class  problem.  Class  one  consists  of  the  KLT  coefficients 
of  sixteen  different  subjects.  Class  two  consists  of  the  KLT  coefficients  of  sixteen  different 
images  of  the  test  subject.  A  Backpropagation  Neural  Network  (BPN)  with  two  layers  and 
two  hidden  nodes  trains  to  differentiate  between  the  two  classes.  To  test  the  classification 
accuracy  of  the  BPN,  each  trmning  vector  is  held  out  once.  The  results,  averaged  over  ten 
trsuls  and  shown  in  Figure  4.14,  demonstrate  the  accuracy  versus  the  number  of  features. 
Five  KLT  coefiUcients  resulted  in  92%  recognition  accuracy. 

Using  the  same  images,  the  single  person  verification  problem  is  run  again  with  a 
Euclidean  distance  classifier.  The  results  shown  in  Figure  4.14  illustrate  the  recognition 
accuracy  versus  the  number  of  features.  Three  KLT  coefficients  resulted  in  93%  recognition 
accuracy. 

Note  that  the  probability  of  error  of  the  Euclidean  distance,  also  called  Nearest 
Neighbor  (NN)  classifier,  approximates  a  Bayesian  classifier  as  the  number  of  test  and 
prototype  vectors  approaches  infinity  (25:83).  Additionally,  Ruck  proves  that  the  BPN 
approximates  a  Bayesian  classifier  which  indicates  that  both  the  BPN  and  NN  classifiers 
are  Bayesian  (17).  The  results  of  Figure  4.14  confirm  this  result. 

Both  tests  indicate  that  for  single  person  verification,  a  small  number  of  coefficients 
provides  relatively  accurate  verification. 

Multi'person  Recognition.  Using  fifty-five  video  taped  subjects,  six  128  by  128 
pixel  frames  of  each  person  are  taken.  One  frame  of  each  subject  is  used  in  the  KLT 
training  set. 

Before  feature  extraction,  the  six  images  are  divided  into  a  test  and  prototype  set. 
The  test  set  consist  of  two  images  which  will  be  compared  to  the  remaining  four  in  the 
prototype  set.  With  the  KLT  coefficients  of  both  the  test  and  prototype  sets  provided  by 
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Figure  4.14.  Results  of  single  person  recognition  versus  number  of  eigenfaces  for  both 
BPN  and  Euclidean  distance  classifiers 


the  routine  recon.c,  a  Euclidean  distance  classifier,  find-min.c,  is  employed  to  compare  the 
test  set  to  the  prototype  set.  A  minimum  distance  implies  a  match. 

The  multi- person  test  first  varies  the  number  of  KLT  coefficients.  Figure  4.15 
demonstrates  the  recognition  accuracy  versus  the  number  of  coefficients.  Ninety-five  per¬ 
cent  recognition  accuracy  is  achieved  with  only  sixteen  coefficients.  Interestingly  enough, 
increasing  the  number  of  coefficients  actually  decreases  accuracy.  This  can  be  readily 
expladned  by  pointing  out  that  the  lower  eigenvalued  eigenvectors  have  less  energy  and 
therefore  these  coefficients  introduce  more  noise  than  information  into  the  classification 

(7). 

A  second  test  is  performed  on  the  same  subjects.  The  test  set-up  is  identical  to 
the  previous  test  except  forty  images  are  used  in  the  KLT  training  set.  Figure  4.16 
demonstrates  the  recognition  accuracy  versus  the  number  of  coefficients  used.  In  this  case, 
ninety-four  percent  accuracy  is  achieved  using  sixteen  coefficients. 

Using  the  Fast  Fourier  Transform  (FFT),  the  same  set  of  Gaussian  windowed  images 
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Figure  4.16.  Results  of  multi-person  recognition  versus  number  of  eigenfaces  for  55  sub¬ 
jects  where  only  40  of  the  subjects  are  included  in  the  KLT  training  set 
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are  transformed,  and  the  FFT  coefficients  are  calculated.  After  Euclidean  distance  classi¬ 
fication,  recognition  accuracy  versus  the  number  of  FFT  coefficients  is  determined.  Figure 
4.17. 


Figure  4.17.  Results  of  multi-person  recognition  versus  number  of  FFT  coefficients  for  55 
subjects 

A  comparison  of  the  KLT  and  FFT  is  shown  in  Figure  4.18.  The  KLT  clearly  out 
performs  the  FFT  for  a  lower  number  of  coefficients.  For  example,  the  FFT  needs  27 
coefficients  to  reach  95%  recognition  accuracy  where  as  the  KLT  only  requires  as  few  as 
16  coefficients.  The  FFT  recognition  performance  approaches  the  performance  of  the  KLT 
as  the  number  of  coefficients  is  increased. 

Finding  Faces  in  Face  Space 

Finding  a  face  in  a  scene  is  the  first  step  in  any  pattern  recognition  system.  The 
concept  of  using  the  KLT  for  segmentation  is  evaluated.  The  first  segmentation  approach 
scans  a  scene.  Each  sub-image  scanned  is  projected  into  face  space.  If  the  Euclidean 
distance  of  the  sub-image  KLT  coefficients  is  less  than  a  certain  threshold,  a  face  is  probably 
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Figure  4.18.  Comparison  of  multi-person  recognition  versus  number  of  coefficients  for 
FFT,  KLT  using  55  training  subjects,  and  KLT  using  40  subjects 


present  and  the  sub-image  is  segmented  from  the  scene.  The  results  of  this  experiment 
indicate  that  faces  can  not  be  found  using  this  technique.  The  technique  is  very  sensitive 
to  head  size  and  lighting.  Worst  of  all,  the  technique  is  computationally  impractical. 

The  second  segmentation  approach  utilizes  the  FFT  of  the  scene  and  of  the  eigenfaces. 
A  correlation  is  performed  beween  the  scene  and  eigenfaces  to  determine  the  location  of 
the  face  in  the  scene.  The  results  of  this  test  indicate  that  faces  can  not  be  found  with 
this  technique.  A  correlation  peak  does  occur  at  the  location  in  the  scene  corresponding 
to  the  face,  but  the  peak  is  not  a  global  extrema  and  many  peaks  have  to  be  tested  before 
a  face  is  confirmed. 

Although  both  of  the  previous  techniques  performed  poorly  because  of  size  and  light¬ 
ing  sensitivity,  the  AFRM  moving  target  segmenter  has  been  shown  to  perform  over  various 
lighting  and  image  size  conditions  (10).  The  AFRM  uses  motion  to  segment  the  image, 
then  the  image  is  normalized  for  brightness  and  size.  With  AFRM  preprocessing,  the 
image  is  properly  normalized  for  a  KLT  based  recognition  system. 
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Another  face  segmentation  technique  demonstrated  by  Tarr,  finds  human  eyes  in  a 
scene  using  a  neural  network  (23).  This  neural  network  based  segmentation  can  also  be 
used  instead  of  moving  target  segmentation.  Tarr  has  also  demonstrated  a  Karhunen-Loeve 
Neural  Network  which  implies  that  a  totally  neural  based  solution  can  implement  many  of 
the  techniques  demonstrated  in  this  thesis  (24). 

Karhunen>Lohve  based  Facisd  Feature  Communicator  for  the  Non-vocal 

The  test  images  are  processed  and  classified  as  in  Figure  3.10.  The  KLT  coefficients 
of  the  test  images  are  compared  to  the  KLT  coefficients  of  previously  processed  prototypes. 
The  nearest  neighbor  classifier  determines  the  class  of  the  image,  i.e.  ‘yes’  or  ‘no’.  The 
system  was  able  to  determine  ‘yes’  tongue  out  or  ‘no’  mouth  open  with  100%  classification 
accuracy. 

The  test  imposed  several  limitations  on  processing.  First,  the  KLT  is  sensitive  to 
size.  So  a  real  system  would  have  to  fix  the  distance  between  the  user  and  camera.  The 
second  limitation  is  segmentation  and  alignment  of  images.  This  limitation  can  be  solved 
in  an  operational  system  by  utilizing  positive  feedback.  Simply  indicate  non-alignment  to 
the  user  and  ignore  all  test  images  that  do  not  fall  within  some  KLT  coefficient  distance 
threshold.  A  non-alignment  can  easily  be  indicated  on  a  computer  screen  and  when  the 
user  is  aligned  properly  the  ‘yes’  or  ‘no’  indications  can  be  processed. 

The  KLT  successfully  determined  ‘yes’  or  ‘no’  in  a  small  test  set  of  four  images.  The 
ability  to  differentiate  between  these  simple  expressions  indicates  that  the  KLT  may  be  used 
as  a  communications  interface  for  the  disabled.  In  other  words,  this  concept  can  be  used 
by  a  disabled  person  to  manipulate  a  computer  menu  system  that  controls  other  functions. 
Similarly,  the  binary  ‘yes’  or  ‘no’  scheme  can  easily  be  coded  into  morse  code  permiting  a 
simple  binary  to  text  conversion.  The  methodology  used  here  should  be  readily  employable 
at  frame  rate  because  the  KLT  coefficients  are  computed  with  only  one  subtraction  and 
two  dot  products  between  the  input  image  and  each  of  the  two  eigenvectors. 
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Karhunen-Lo^ve  Axis  System 

The  purpose  of  this  test  is  to  demonstrate  correlation  in  Karhunen-Loeve  space  and 
to  develop  an  axis  system  or  technique  to  align  laser  scanner  head  images.  The  first 
eigenvector,  Figure  4.19,  is  the  basis  vector  with  the  greatest  contrast  or  variance  among  all 
the  images.  Therefore,  an  eigen  based  axis  system  takes  into  account  the  largest  variances 
in  all  of  the  training  set  population. 


Figure  4.19.  An  example  of  an  eigenvector  resulting  from  a  training  set  consisting  of  head 
scanner  data. 


The  graph  of  the  projections  is  shown  in  Figure  4.20.  The  projections  demonstrate 
that  as  the  image  approaches  the  origin  of  the  axis  system,  the  projection  value  approaches 
a  peak.  The  peak  corresponding  to  the  axis  center  is  global.  After  the  extrema  is  found, 
the  image  is  shifted  at  one  pixel  increments  to  find  the  precise  peak.  The  projections  at 
one  pixel  resolution  about  the  axis  origin  are  shown  in  Figure  4.21.  The  peak  at  (5,-32) 
corresponds  to  a  shift  of  the  test  image  5  to  the  right  and  32  down. 

To  compare  the  result  of  the  Karhunen-Loeve  correlation  with  the  FFT  based  cor¬ 
relation,  the  same  image  is  correlated  with  the  first  eigen  image  which  was  used  in  the 
previous  test.  The  FFT  based  correlation  indicated  a  correlation  peak  at  (5,-32).  This  is 
in  agreement  with  the  results  obtained  by  the  previous  Karhunen-Loeve  approach. 
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Figure  4.20.  The  top  graph  shows  the  projection  of  the  shifted  image  into  Karhunen- 
Loeve  space  using  only  the  first  eigenvector.  The  bottom  graph  details  the 
peak  location  which  indicates  alignment  between  the  shifted  test  image  and 
eigenvector. 
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Figure  4.21.  Projection  of  the  shifted  test  image  with  one  pixel  resolution  into  Karhunen- 
Loeve  space  using  only  the  first  eigenvector.  The  minimum  represents  the 
correlation  peak  indicating  alignment  between  test  image  and  eigen  image. 


The  results  indicate  that  correlation  can  be  performed  with  the  Karhunen-Loeve 
coefficients,  and  the  coefficients  can  be  used  to  align  images.  The  similarity  in  the  location 
of  the  correlation  peak  with  Fourier  based  correlation  lends  confidence  to  the  result.  Recall 
that  the  first  eigenvector  represents  an  image  with  the  greatest  amount  of  variance  or 
information  in  its  pixels.  The  eigenvector  provides  the  best  image  to  correlate  against 
because  the  eigenvector  represents  the  largest  variance  or  worst  case  change.  In  short, 
the  first  eigenvector  provides  the  best  correlation  with  the  largest  population.  This  is  the 
justification  for  the  use  of  the  eigen  images  as  a  basis  for  an  axis  system.  Further  evaluation 
needs  to  determine  the  robustness  of  an  eigen  based  axis  system  over  a  larger  population 
of  test  images.  The  use  of  the  eigen  images  with  either  an  FFT  or  KLT  based  correlator 
seems  to  be  equivalent. 

The  significant  results  of  this  research  are  the  demonstration  of  correlation  in  Karhunen- 
Loeve  space  and  its  equivalence  to  traditional  FFT  correlation.  Second,  the  use  of  eigen 
images  as  a  basis  for  an  axis  system  represents  a  logical  choice  over  traditional  feature  based 
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techniques,  e.g.  an  axis  system  based  on  the  planes  traversing  the  nose  and  ears.  Further 
research  will  evaluate  the  results  over  a  larger  data  set  and  evaluate  the  use  of  a  conjugate 
gradient  approach  to  Karhunen-Loeve  correlation  techniques  which  would  eliminate  the 
need  for  as  many  shifts  to  find  the  global  extrema. 

Conclusion 

The  most  significant  result  of  this  chapter  is  the  demonstration  of  the  recognition 
superiority  of  the  KLT  over  the  FFT  or  the  AFRM  feature  based  approach.  With  only 
16  coefficients,  the  KLT  provided  95%  face  recognition  accuracy  for  a  population  of  55 
subjects.  The  40  subject  test  demonstrated  KLT  generalization  for  recognition.  This  is 
significant  because  the  poor  reconstruction  of  non-training  subjects.  Figure  4.13  ,  indicates 
that  for  reconstruction  purposes  the  KLT  did  not  generalize.  In  short,  this  demonstrates 
that  the  KLT  is  good  for  face  recognition  but  not  sufficiently  adequate  for  reconstruction. 

This  chapter  first  presented  the  results  of  code  testing  and  then  presented  the  results 
of  the  thesis.  The  chapter  then  presented  the  reconstruction  results  for  both  large  and  small 
population  training  sets.  Thirdly,  the  chapter  provided  results  of  face  recognition  using 
the  KLT  and  FFT.  The  last  sections  presented  the  results  the  KLT  based  face  finder,  the 
facial  feature  communicator,  and  the  Karhunen-Loeve  axis  system.  The  significant  results 
and  conclusion  are  presented  in  the  next  chapter. 
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V.  Conclusions 


Introduction 

The  purpose  of  this  thesis  is  to  find  good  features  for  machine  recognition  of  human 
faces.  This  thesis  investigated  the  Karhunen-Lohve  Transform  (KLT)  and  demonstrated 
its  performance.  A  summary  of  the  significant  results  follows. 

Summary  of  Significant  Results 

This  thesis  developed  efficient  KLT  algorithms  in  ANSI  C.  The  algorithms,  which  are 
documented  in  Appendix  A,  implemented  a  Karhunen-Loeve  approximation  that  simplified 
the  solution  to  a  simultaneous  equation  from  degree  to  M.  For  example,  for  128  by 
128  pixel  images  in  a  KLT  trmning  set  of  40  images,  =  16384  and  M  =  40.  This 
approximation  which  is  only  valid  if  >>  M  greatly  simplifies  and  hence  speeds  the 
calculation  of  the  KLT  (27). 

This  thesis  demonstrated  the  image  coding  and  compression  capabilities  of  the  KLT. 
Large  and  small  populations  of  subjects  were  used  in  the  KLT  trmning  set.  For  a  small 
population  of  six  subjects,  the  reconstruction  quality  of  the  KLT  encoded  image  was  as 
good  as  the  original  image.  The  KLT  was  able  to  compress  a  256  by  256  pixel  image  into 
three  coefficients.  This  represents  a  compression  of  approximately  21,000:1.  For  a  larger 
population  of  55  subjects,  the  reconstruction  of  the  128  by  128  encoded  images  was  reason¬ 
ably  good  with  only  twenty  coefficients.  This  represents  a  compression  of  approximately 
800:1.  Both  of  these  dramatic  compressions  come  at  the  expense  of  a  constrained  subject 
population  and  pre  and  post  processing  requirements. 

The  most  significant  achievement  of  this  thesis  is  the  dramatic  face  recognition  ac¬ 
curacies  obtcuned  using  the  KLT  coefficients.  The  KLT  coefficients  provided  93%  single 
person  verification  with  as  little  as  three  coefficients.  For  multi-person  recognition,  the 
KLT  coefficients  achieved  95%  accuracy  with  as  little  as  16  coefficients  and  a  population 
of  55  subjects.  In  comparison,  the  FFT  with  the  same  test  data  achieved  85%  with  16 
coefficients.  These  results  are  markedly  improved  over  the  inconsistent  AFRM  accuracy 
which  varied  from  56%  to  90%  for  a  population  of  50  subjects  and  15  coefficients.  The 
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KLT  outperforms  the  FFT  and  the  AFRM.  The  KLT  is  sensitive  to  changes  in  head  size 
and  position,  but  the  FFT  based  techniques  also  suffer  from  this  problem. 

This  thesis  demonstrated  that  a  transform  may  be  good  for  recognition  but  not 
adequate  for  reconstruction.  The  Karhunen-Loeve  Transform  demonstrated  very  good 
recognition  accuracy  for  subjects  out  of  the  KLT  training  set.  This  is  significant  since  the 
KLT  provided  poor  reconstructions  of  out  of  training  set  subjects. 

Another  significant  result  of  this  effort  is  that  the  software  written  for  faces  can  also 
be  easily  used  for  any  other  image  recognition  application.  The  software  written  for  this 
thesis  was  utilized  by  Singstock  to  recognize  tanks,  trucks,  jeeps,  and  towers  in  forward 
looking  infrared  imagery  (20).  Singstock  developed  eigentanks,  eigentrucks,  eigenjeeps, 
and  eigentowers  as  a  basis  set.  The  eigen  images  were  successful  in  classifying  tank,  truck, 
jeep,  and  tower.  Figure  5.1  displays  a  Singstock  eigentank. 


Figure  5.1.  The  Singstock  eigentank  used  to  recognize  tanks. 


A  related  application  to  face  recognition  is  facial  expression  recognition.  A  test 
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was  conducted  to  differentiate  facial  expressions.  The  motivation  behind  the  test  was  to 
evaluate  the  possibility  of  developing  a  facial  feature  communications  interface  for  a  child 
with  cerebral  palsy.  The  child  is  severely  disabled  and  can  not  speak.  The  child  does  have 
control  of  facial  muscles,  therefore  differentation  between  tongue  out  and  mouth  open 
provides  a  means  of  communication.  The  Karhui.^:?-Loeve  system  was  able  to  determine 
‘yes’  tongue  out  or  ‘no’  mouth  open  with  100%  classification  accuiury. 

The  last  related  application  is  in  the  anthropometric  community.  The  anthropometric 
community  needs  an  axis  system  for  head  laser  scan  data.  The  head  laser  scan  data  is 
basically  an  unwrapped  picture  of  the  head.  Before  anthropometric  data  can  be  extracted, 
the  image  must  be  centered  to  simplify  the  measurement  process.  This  thesis  evaluated 
the  use  of  a  Karhunen-Loeve  based  axis  system.  The  results  indicate  accurate  placement 
of  laser  scanner  heads  with  the  KLT  axis  system. 

Conclusions 

In  conclusion,  this  thesis  demonstrated  superior  face  recognition  using  the  KLT.  The 
KLT  approximation  implemented  in  this  thesis  has  removed  the  computational  burden 
associated  with  Karhunen-Loeve.  Also,  note  that  the  results  using  a  simple  first  nearest 
ndghbor  classifier  were  95%.  To  attain  even  greater  accuracy  a  more  elaborate  classification 
scheme  such  as  a  neural  network  could  be  implemented.  Additionally,  a  multiple  look  or 
multiple  frame  collection  of  each  subject  would  also  improve  the  recognition  accuracy  above 
95%. 

The  primary  limitation  of  the  Karhunen-Loeve  is  scale  and  position  sensitivity.  There 
are  two  solutions  to  this  problem.  First,  the  AFRM  segmenter  sucessfuUy  provides  scale, 
contrast,  and  lighting  normalized  images.  Therefore,  a  Karhunen-Loeve  based  recognition 
system  could  utilize  the  AFRM  segmenter.  A  second  approach  is  to  transform  the  image 
into  a  scale,  position,  and  rotation  invariant  space  as  demonstrated  by  Kobel  and  Martin 
(8). 

Lastly,  the  software  written  for  this  thesis  is  flexible  and  the  techniques  developed 
here  have  been  applied  to  problems  associated  with  tank  recognition,  anthropometry,  and 
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communications  for  the  severely  disabled. 


Future  Areas  of  Research  and  Recommendation 

The  prime  recommendation  for  future  work  would  be  to  integrate  the  work  done  on 
this  thesis  with  the  AFRM.  This  would  take  advantage  of  the  AFRM  camera  and  segmenter 
and  utilize  the  Karhunen-Loeve  techniques  developed  here.  Additionally,  a  fruitful  area  of 
research  would  be  to  implement  a  total  neural  network  approach  to  the  face  recognition 
problem.  Finally,  the  human  face  has  a  unique  color  signature  and  future  research  should 
investigate  the  use  of  color  cameras  for  segmentation  and  recognition. 

General  Summary 

This  chapter  provided  a  summary  of  the  significant  results  as  well  as  recommen¬ 
dations  for  future  research.  Chapter  1  provided  background  and  defined  the  problem, 
objectives,  questions,  methodology,  standards,  and  scope  of  the  research.  Chapter  2  pro¬ 
vided  auditional  background  and  a  review  of  current  research.  The  methodology  used  to 
solve  the  problem  was  presented  in  Chapter  3.  The  results  of  the  research  were  presented 
in  Chapter  4.  Appendix  A  lists  all  the  code  utilized  for  this  thesis. 
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Appendix  A.  Thesis  Source  Code 


This  appendix  contains  a  listing  of  all  code  developed  and  or  used  for  this  thesis. 
This  code  is  presented  as  is,  and  no  claims  are  made  as  to  suitability  for  other  applications. 
Portions  of  this  code  are  copyrighted  by  the  Air  Force  Institute  of  Technology  and  by  Pedro 
F.  Suarez.  Such  portions  are  indicated  as  such. 


kl_transform2  .c 

NAME:  kLtransforml.c 
INVOKED: 

DATE:  24  June  1991 
DESCRIPTION: 

SUBROUTINES  CALLED: 

FUTURE  MODIFICATIONS/BUGS: 

#include  <stdio.h> 

#include  <math.h> 

^include  <string.h> 

#define  SQ(A)  (A*A) 

main(argc,argv) 
int  argc; 
char  «argv[]; 

{ 

FILE  -strain,  4>facein,  ■i-fout; 

FILE  *face^vg,  ^tempfUe,  +fevex,  ♦feval; 
int  ij,  N,  k,  M,  nrot,  atoi(); 

float  *«matrix(),  ♦vector(),  ♦*A,  4<4>A.trans,  ♦♦u,  ♦♦L,  **v,  ^laverageJace; 
float  *d,  temp,  *mag; 

void  free_vector(),  free_matrix(),  eigsrt(),  jacobi(),  mat.col_mag(); 
char  fllename[81],  h^; 


if  (argc  /  4)  { 

printf("!!!  The  comnand  line  should  be  !!!:\n\n  kl_tTansfoTm2  trainingfile 
size  i_in_train.set\n\n"); 
exit(O); 

} 

/**:t*\  Set  Up  Files 

if  ((train  =  fopen(argv[l],  "r"))  ==  NULL) 


A-1 


{ 

prmtf("I  can't  open  the  training  file"); 
exit(— 1); 

} 

/**** 

printff”///  luput  total  image  SIZE(N):”); 

scaB{(”%d”,  &N); 

pnntf("\n"); 

pnntf(”!!!  Input  the  number  of  training  faces  (M):”); 

scanf(”%d”,  &M); 

printf(”\n”); 

****! 

N=atoi(argv[2]); 

M=atoi(argv[3]); 


/*  dynamically  allocate  memory  */ 

A.trans  =  matrix(l,M,l,N); 

A  =  matrix(l,N,l,M); 

average  Jace  =  vector  (1,N); 

L  =  matrix(l,M,l,M); 
d  =s  vector(l,M); 

V  =  matrix(l,M,l,M); 
mag  =  vector(l,  M); 

/****  initalize  matrix  and  vectors  ****/ 
for(j=lu<Mu++) 
for(i=l;i<Nd++)  { 

A_tran8[j](i] = A[i]|j] =average  Jacep] =0.0; 

} 

for(k=l;  k<M;  k++){ 

f8canf(triun,  "XsNn",  fUename); 
/***printf(”%s\n”f  filename);  ***/ 

facein  =  fopen(iUename,  "r"); 


for(j=l;j<Na++){ 

f8canf(facein,"Xf\n",&A[j][k]); 

} 

fdo8e(facein); 

/***  nrint/f  "first  value  of  file  %d  =  %f\n”,  It,  A[l}[k]);  ♦♦*/ 
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!***  Normalizing  Data  by  dividing  by  255  ***/ 


/****  for(j=lu<Mu+i-) 
for(i=l;i<N;i++)  { 
A[iJ[iJ=A[i][j]/255; 

} 

*******! 


/******rif***0CaIcttlate  Average  Face 

printf('M  !!  CALCULATE  AVERAGE  FACE\n  "); 
face^vg=fopeii("avg_f  ace .  dat" 
for(i=l;i<N;i++){ 
temp=0.0; 
for(j=lu<My++) 
temp=temp+A(i][j]; 
averageJace[i]=temp/M; 
fprintf(face Jivg,"Xf  \n" ,  average_face[i]); 

} 

fclose(face^vg); 


/***  printf(”!!!  SUBTRACTING  OFF  AVERAGE  FACE  \n”);  ***/ 
for(j=la<M;j++) 
for(i=l;i<N;i++){ 

A  [i]  [j] = A  [i]  p]  -  aver  age  Jace[i] ; 

} 


/**♦♦  CREATING  A  TRANSPOSE  ***/ 
printf("!!!  CREATING  TRANSPOSE  MATIX  W); 

for(j=ly<My++) 

fbr(i=l;i<N;i++){ 

A-tran8p](i]=A[i]p]; 

} 


/***  font =fopen (”l.dat”,”w”);**^ 

printf("!!!  Multiplying  A  trans  and  A  to  get  L:\n"); 
for(i=l;i<M;i++) 
for(j=lu<M;j++)  { 
temp=0.0; 
for(k=l;k<N;k++){ 


A-3 


temp=temp+A-trans[i][k]*A[k][j]; 

/**fprintf(fout,”%f\n”,  temp);**/ 
L[il[j]=temp; 

/**  pnat{(”!!!  Writing  Output  \a”);  ***/ 


/****  priat{(”!!!  L  Matrix  WRITTEN  TO  12.DAT  \b”); 

fclose( lout);****/ 


/***  (”!!!  FREE  MATRIX  A.TRANS  \b”);  ***/ 
free_matrix(A-tran8, 1,  M,  1,  N); 


printfC'!!!  doing  jacobian  of  L  \n''); 

ja«obi(L,M,d,v,  &nrot); 

/*** print f(”%d\n”,  nrot);***/ 

/***  printf(”!!!  doing  eigsrt  of  v  \n”);**^ 
eig8rt(d,v,M); 


printfC'  *  *  *  Writing  aigenvvalues  \n"); 
feval=fbpen("aigan.val"  ,"b"  ); 
for(j=ly<Mij++)  { 

fprintf(feval,"Xf\n",  dLj]); 
fdo8€(feval); 


printf("!!!  Writing  aiganvectors  \n"); 

fevex=fopen("aigan_vac","w"); 
for(i=l;i<M;i++){ 
for(j=ly<My++)  { 

fprmtf(fevex,"Xf  \n",  v|j][i]); 


fdo8e(fevex); 


/**  print  f(”!!!  Initailizing  Eigen  face  Matrix  \n”);  **/ 
u  =  matrix(l,N,l,M); 
for(k=l;k<M;  k++) 
for(j=ly<N;  j++) 
u[j][k]=0.0; 
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printf(''!M  Calcualting  aiganfac*  \n"); 


for(k=l;k<M;  k++){ 
for(i=l;i<M;  i++){ 
for(j=l*J<N;  j++) 

«DlM=vli][k]*AD]ri]4-uLj][k]; 

} 


/***iindmg  magnitude  of  eigenface  ***/ 
/**  matjcoI-mag(u,  N,  M,  mag);  ***/ 


printf("!!!  Opening  train. out  file  for  Eigenfaces  \n"); 
tempfUe  =  fopen( "train. out",  "u"); 


h=48; 

1=48; 

strcpy(fUename,  "eigenface"); 

for(k=l;  k<M;  k++){ 

if(l  57)  1++; 
el8e{ 

1=48; 

h++; 

} 


fprintf(tempfUe,"  dename); 
fprintf(tempfUe,"XcXc"  ,h4); 
fprmtf(tempfUe,"X8\n","  .dat"); 

} 

fclose(tempfUe); 


printf("!!!  Writing  eigenface  \n"); 
tempfUe  =  fopen(" train. out",  "r"); 
for(k=l;  k<M;  k++){ 

fscanf( tempfUe,  "XaVn",  filename); 
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facdn  =  fopeii(fileiiajne, 

/***  pnBtf(”%8\n”,  aieBame);  ***/ 

/***  fpriBt{(faceiB,  mag[k]);  ***f 

for(j=la<Ny++){ 
ft(rintf(facein,"Xf  \n'',  u[i][k]); 

} 

fclo8e(facein); 

} 


fclo8e(tempfile); 

l***priBt^”VA  FREEING  A  MATRIX 
freejnatrix(A,l,N,l,M); 

free-matrix(u,l,N,l,M); 

} 

void  matjcoI^ag(u,  N,  M,  mag) 
float  4>*u,  magQ; 
int  N,M; 

{float  b; 
int  k,  j; 

double  sqrt  (); 

for(k=l;  k<M;  k++)  { 
b=0; 

for(j=l;  j<N;  j++) 
b=u(j][k]  *  u[j][k]  +  b; 
mag[k]  =  8qrt(  (double)  b); 

} 

} 


recon.c 

NAME:  recon.c 
INVOKED:  reconstructs  faces 
DATE:  24  July  1991 
DESCRIPTION: 

SUBROUTINES  CALLED: 

FUTURE  MODIFICATIONS/BUGS: 

/***  ^mclude  <device.b>  ♦♦/ 

#  include  <stdio.h> 

#indude  <math.h> 

#include  <string.h> 
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#define  SQ(A)  (A*A) 
maiii(argc,argv) 
int  argc; 
char  ^argvfl; 

{ 

FILE  4>facel,  ^eigenin,  *fout,  *fweights,  *fspectnun,  strain; 

FILE  *face^vg; 
int  ij,  N,  k,  M,  atoi(); 

float  *vector(),  4>*matrix(),  ^averagejace,  max; 

float  *d,  temp,  *pedro,  *reconface,  *w,  *1,  temp-mag,  mag; 

void  free_vector(); 

void  rescaleO; 

void  vecjnag(); 

double  sqrtQ,  fabsO; 

char  filename[81],  *8trcpy(),  infllename[81],  outfllename[81],  ext[10],  junk[81],  ngfile[81]; 
if(argc  7^  5){ 

printf("!!!  The  command  line  should  be  !!!:\n\n  recon  infile  ngfile  imagesize 
numbarracon\n\n" ); 
exit(O); 

} 

Set  Up  Files 

8trcpy(ext,  ".gra"); 

8trcpy(infilename,argv[l]); 
strcat(infilename,  ext); 

8trcpy(ext,  ".rec"); 
strcpy(outfllename,  argv[l]); 

8trcat(outftlename,  ext); 


if  ((facel=fopen(infilename,"r"))  ==  NULL){ 
printf("I  can't  open  the  input  file"); 
exit(-l); 

} 


if  ((fout=fopen(outfilename,"u"))  ==  NULL){ 
printf("I  can't  open  the  output  file"); 
exit(-l); 

} 

N  =  atoi(argv[3]); 

M  =  atoi(argv[4]); 


A-7 


printi(”!!!  Input  total  image  SIZE(N):”); 

aamf(”%d”,  &N); 

print 

printf(”!!!  Input  the  number  of  training  faces  to  reconstruct  from  faces  (M):”); 

8canfC%d”,  &M); 

printf(”\n”); 


/*  dynamically  allocate  memory  */ 

u  =  matrix(l,N,  1,M); 
pedro  =  vector(l,  N); 
average-face  =  vector(l,  N); 
reconface  =  vector(l,  N); 
w  =  vector(l,  M); 

I  =  vector(l,  N); 

!****  initalize  and  vectors  ****/ 
for(j=l;  j<M;  j++) 
for(i=l;i<N;i++) 

w[jj=u(i]  [j] =I[i] =pedro[i] =reconface[i]=average_face[i] =0.0; 


/*«*«  printf(”!!!  Opening  and  Reading  Pace  to  be  Reconstructed:\n  ”);  *♦*+/ 

for(j=lij<Ny++) 
f8canf(face  1 , "  Xf  \n" ,  &pedro[j] ) ; 
fclose(facel); 


/*  printf(”!!!  Opening  and  Reading  Face  Average  Face  (avgJace.dat):\n  ”);*/ 

face^vg=fopen("avg_lace .dat",  "r"); 
for(j=lu<N'J++) 

f8canf(facejivg,  "XfVn",  &averageJace[j]); 
fclo8e(face-avg) ; 

/**  printf(”!!!  Using  %d  eigenfaces  from  file  train.out;\n  ”,  M);  **/ 

train  =  fopen(  "train. out",  "r"); 

for(j=l;  j<M;  j++){ 

f8canf( train,  "XaXn",  filename); 

/**  printf(”%s\n”,  filename), *  **/ 
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eigenin  =  fope]i(iUename,  "r"); 


for(i=l;i<N;i++){ 
fscanf(eigemn Xf  \n"  ,&u[i]  [j] ) ; 

} 

fdose(eigenin); 

/»♦  printf(”first  value  of  file  %d  =  %f\n”,j,  u[l}[j]);  **/ 


fclo8e(train); 


printfC'!!!  SUBTRACTING  OFF  AVERAGE  FACE  \n''); 
for(i=l;i<N;i++){ 

I[i]=  pedro[i]  —  average  Jace[i]; 

} 


/*  printf(”!!!  Calculating  Weights  based  on  %d  eigen faces\n”,  M);  ♦*/ 

for(j=l;  j<M;  j++) 
for(i=l;  i<N;  i++){ 

wp]  =  u[i]p]*  I[i]+  wp]; 

} 


fspectrum  =  fopen("8pectr\m.out",  "a"); 

fweights  =  fopen(argv[2],  "a"); 
for(i=l;  i<M;  i++){ 

fprintf(fweights,  "Xf  ",  w[i]); 
fprintf(f8pectrum,  "Xd  Xf\n",  i,  w[i]); 

} 

fprintf( fweights,  "XsXn",  argv[l]); 

fclo8e(fweight8); 
fclo8e( fspectrum); 

/***  Normalize  weights  by  dividing  by  absolute  max  *****/ 
max  =  0; 

for(i=l;  i<M;  i++){ 

/***  printf(”%f\n”,w[i]);**^ 
if(fab8(  (double)  w[i])  >max){ 
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max=  fabs(  (doable)  w[i]); 

} 

} 

for(i=l;  i<M;  i++) 
w[i]=w[i]/inax; 


/♦♦*♦  FACE  RECONSTRUCTION  ****/ 


printf("!!!  Reconstructing  fron  XM  facesW,  M); 

for(j=l;  j<M;  j++) 
for(i=l;  i<N;  i++) 

reconface[i]  =  w[j]  *  u[i][j]+  reconfacefi]; 


printf("!!!  Adding  mean  back  on  Reconstructing  facesXn"); 
for(i=l;  i<N;  i4-+) 

reconfacefi]  =  reconfacefi]  +  (averageJacep]); 
re8caIe(reconface,  N); 

printf("!!!  Writing  Reconstructed  Pedro  \n  ”); 

for(j=lJ<Nij++){ 

fprintf(fout,  "X4.0f\n",  reconface[j]); 

/♦♦♦  fwrite(&TecoBface[j],  sizeof(reconface[j]),l,fout); 

fclo8e(fout); 


} 


void  re8cale(output,  N) 
float  output  []; 
int  N; 

{ 

int  NEW_MAX,  NEW3iIN,  i,  j,  count; 
float  min,  max; 

NEW3IAX  =  255; 

NEW_MIN  =0; 
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/***  priat{(”\n!!!  FINDING  MAX  n!\a\B”);  ***/ 


/**  Check  for  the  max  and  min  value  in  the  data  **/ 

min=output[l]; 

max=outpat[l]; 

cottiit=0; 


for(j=l;  j<N;  j++){ 

if(oatput  [j]  >max)  { 

/**  printf(”%f\a”,  max);***/ 
max=output[j]; 

} 

if(output[j]<min){ 

min=output[j]; 

/♦*  priBtf(”%f\n”,  min);  **/ 

} 

count++; 

} 

/♦*  pnatf(”\n  SAMPLES  =  %d\n”,  count); 
printf(”  max  =  %f  min  =  %f\n”,  max,  min);  **/ 


!***  Now  translate  data  and  write  to  output  hie  *♦<•/ 
for(i=l;  i<N;  i++){ 

output[i]  =  ((output[i]-min)*(NEW_MAX-NEW_MIN)/(max-inin)  +  NEW31IN); 

} 

/♦♦  printf(”!!!  All  Done  !!!\n\n”);  **/ 

} 

void  vecjnag(input,  N,  mag) 
float  mputQ,  ^‘mag; 
int  N; 

{float  b; 
int  i; 

double  sqrtO; 
b  =  0; 


for(i=l;  i<N;  i+4-){ 

b  =  input[i]  *  input[i]  +  b; 


} 


All 


b  =  sqrt(  (double)  b); 

/***  printf(”mag  for  u  =  b);  ♦♦♦</ 

«mag  =  b; 

} 


find-min.c 

\*  This  program  implements  a  First  Nearest  Neighbor  Network. 

Date:  Aug  1991 

Written  By:  Pedro  F.  Suarez 

Invoked:  findjuin  trainfUe  testiile  #-train  #-test  ^.features 

where  trsuniile  is  a  file  containing  the  prototype  vectors,  testiile  is  a  file  containing  the  test 
vectors,  ^.trsun  number  of  training  vectors  in  trsdnfile,  number  of  test  vectors  in 

testiile,  #  Jeatures  in  each  vector «/ 

#include  <device.h> 

#include  <stdio.h> 

^include  <math.h> 

#include  <string.h> 

#  define  number  blocks  300 
#define  knumber  4 

#define  SQ(A)  (A*A) 
main(argc,argv) 
int  argc; 
char  <»argv[]; 

{ 

FILE  ♦ftest,  '^ftrain,  ♦fout; 

float  distance,  <«>average,  «sd,  temp,  4<vector(),  *<t>matrlx(),  number  .correct; 
int  number-train,  number  Jeatures,  junk,  temp2,  number  .test,  j,  i,  min,  k,  atoi(); 
double  sqrt(); 

struct  block 

float  ♦♦feature.vecl;  /♦  container  for  feature  vectors  ♦/ 
float  «mean.vec; 

char  name[20],'  /*  container  for  name  */ 

float  distance;  /»  contain  distance  between  test  and  feature  *f 

}; 


struct  block  atest[numberblocks];  /♦  reserve  test  blocks  ♦/ 
struct  block  btrain[numberblocks];  /♦  reserve  train  blocks  ♦/ 

if(argc  6  ){ 

printf("!!!  The  conmand  line  should  be  !!!:\n"); 
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printf("find_min  trainfile  testfile  i.train  f.tast  i_features\n"); 
exit(O); 

} 

printf("  open  files  \n"); 


if  ((ftrain=fopen(argv[l],"r''))  ==  NULL){ 
printf("I  can't  open  the  train  file"); 
exit(-l); 


if  ((fte8t=fopen(argv[2],"r"))  ==  NULL){ 
printf("I  can't  open  the  test  file"); 
exit(-l); 


printf("files  opened  \n"); 
number_train=atoi(argv[3]); 
number  Jeature8=atoi(argv[5]); 
number  . test  =  atoi(argv[4]); 


printf("  initialize  b  block  \n"); 
temp2  =  number.irain/knumber;  /* 
for(i=l;  i<temp2;  i++){ 

btrain/j/.feature.vec  =  matrix(l,kaumber,  1,  number Jeatures); 
btrainfij.mean.vec  =  vector(l,  number Jeatures); 
btraittfij.distance  =  0.0; 

} 

print f(”  initialize  a  block  \n”); 

for(i=l;  Knumber.test;  i++){ 

atest[i].feature.vec=vector(l,  number  Jeatures); 
ate8t[i]. distance  =  0,0; 

} 

printff”  create  vector  \n”); 

average=vector(l,  number  Jeatures); 
sd=vector(l,  number  Jeatures); 

print f(”  init  vector  \n”); 

for(i=l;i<  number  Jeatures;  i++ ) 
average[i ]=sd[i ]=0. 0; 


printf(”  read  in  test  vector  \n”); 
for(i=l;i<numb€rJest;i-f--i-){ 
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for(j=:l;  j<aumberJeature8ij++)  { 
fscan{({test,  ”%f”,  &temp); 
atest[i].feature.vecjj]=temp; 

} 

tscanf({test,  ”%s\a”,  atest[i].name); 

printf(”  read  in  train  vector  \n”); 

for(i=l;  i<temp2;  i++) 

{ 

for(j=l;  j<numberJeatvres‘J++ ) 

{ 

for(k=l;  k<knumber;  k++) 

{ 

tscanf(ftrajn,  iitemp); 
btrain[i].feature.vec[j,k]=temp; 

} 

} 

/*  Note  that  the  same  name  is  used  for  all  k  training  vectors  per  class.  I  did  this  because 
we’re  only  interested  in  making  sure  the  right  class  wins,  and  my  pathetic  C  skills  ran  out 
at  this  point. 

f8canf(ftrajn,  "XsVn",  btrain[i].name); 

} 

printf("  statistically  normalizs  \n''); 

for(j=l;  j<numberJeature8;  j++) 
for(i=l;  i<number_te8t;  i++) 

average[j]+=ate8t[i].feature_vec[j]/(number.test  +  number.train); 

for(j=l;  j<temp2;  j++) 

for(i=l;  i< number.train;  i++) 
for(k=l;  k<knumber;  k++) 

{ 

average[j]+=btrain[i].feature_vecp,k]/(numberJ;est  +  number-train); 

} 

for(j=l;  j<numberJeature8y++) 
for  (i= 1  ;i  <  number  .te8t  ;i+ + )  { 

sdjj]  +=  (ate8t(i].feature_vec[jj  -  average[j])*(atest[i].feature_vec[j]  -average[j]); 

} 

for(j=l;  j<temp2J++) 

for(i = 1  ;i  <  number.train  ;i+ + ) 
for(k=l;  k<knumber;  k++) 

8d[j]  +=  (btrain[i].feature-vec(j,k]  -  average[j])*(btrain[i].feature-vec[j,k]  -  average[j]); 
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} 


for(j= 1  y  <number  Jeaturesy ++ ){ 
sdp]  =  8qrt(  (double)  sd^]); 

sdp]  =  ((double)  l/(number_test +number .train— 1))  ♦  sdQ]; 

} 

fout=fopen(''data_norm.out","w"); 
for(j=l;  j<numberJeature8y++) 
fprintf(fout,"Xf  '',averageij]); 

fprintf(fout,"avarag«\n  "); 

for(j=l;  j<number  Jeature8y++) 
fprintf(fout,"Xf  ",8d[j]); 

fprintf(fout,"8d\n  "); 

for(i=l;  i< number-test;  i++){ 
for(j=l;  j<number  Jeaturesy++){ 
if(8d[j]  ^  0) 

atest[i].feature-vec[j]=  (atest[i].feature_vec[j]  -  average[j])/sd[j]; 
else  ate8t[i].feature-vec(j]=  (atest[i].feature-vec[j]  -  average^]); 

fprintf(fout,"%f  "  ,atest[i].feature.vec[j]); 

fprintf(fout Xs \n  "  ,atest [i]  .name) ; 

} 

for(i=l;  i<number_train;  i++){ 
for(j=l;  j<number  Jeaturesy+-i-){ 

if(8d[j]^  0) 

btrain[i].feature_vec[j]=  (btrain[i].feature.vec[j]  -  average(j])/sd[j]; 
else  btrain[i].feature_vec(j]=  (btrjdn[i].feature_vec[j]  —  average[j]); 

fprintf(fout,"Xf  ",btrain[i].feature-vec[j)); 

fprintf(fout,"X8\n  ",btrain[i].name); 

} 


fclose(fout); 

number-corTect=0; 

for  ( k= 1  ;k  <  number  .test  ;k ++)  { 

for(i=l;  i<number-train;  i++){ 
temp=0; 
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for(j=l;  j<numberJ'eaturesy++)  { 

temp  =  (atest[k].feature_vecp]  -  btrain[i].feature_vec[j])  *  (atest[k].feature_vec[j] 
btrain[i].feature_vec[j])  +  temp; 

btrsun[i]. distance  =  sqrt  (  (double)  temp); 


temp  =  btrun[l].distance; 
min  =1; 


for(i=l;  i<number .train;  i++){ 
if(btrain[i] .distance  <  temp){ 
temp  =  btrmn[i].distance; 
min  =  i; 

} 


} 

printf("the  matching  face  is  for  Xs  is  Xs  with  distance  of  Xf\n", 
atest[k]  .name,btrain[min]  .name,btrain[min]  .distance); 

} 

} 


gwind.c 

/♦♦♦*♦♦♦♦**♦♦♦♦*♦♦♦*♦♦♦**♦♦♦♦♦*♦♦♦♦*******♦♦♦♦*♦♦♦♦*♦♦*******♦♦♦ 
This  routine  takes  an  image  by  a  guassian  window. 


written  by:  Pedro  F.  Suarez 
29  July  91 


#iuclude  <device.h> 
#include  <stdio.h> 
#include  <math.h> 
#deiine  pi  3.1416 
#define  SQ(a)  a«a 

main(argc,argv) 
int  argc; 
char  ♦argvQ; 

{ 
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FILE  *fin,  *10111,  *fo; 

float  **p,  **w,  **matrix(),te8t,  xmean,  xvar,  ymean,  yvar,  normal; 

int  row,  col; 

int  i,  j,  k,  count; 

int  outval,  tempJ,  tempj; 

float  inval,  max,  min; 

double  exp(); 

void  rescaleO; 

if(argc  3){ 

printf("!!!  The  conmand  line  should  be  !!!:\n\n  gauseind  infile 
exit(O); 

} 


if  ((fin=fopen(argv[l],"r"))  ==  NULL){ 

printf("I  can't  open  the  input  file"); 
exit(-l); 

} 

if  ((fout=fopen(argv[2],"w''))  ==  NULL){ 

printf("I  can't  open  the  output  file"); 
exit(-l); 

} 

row  =  col  =128; 
xmean  =  65; 
xvar  =  25; 
ymean  =  55; 
yvar  =  30; 

/****** 

priat{(”Input  the  number  of  rows  and  cols  in  Row  Col 
scanf(”%d  %d”,  krow,  &col); 

priBtf(”%d  %d\n”,  row,  col); 

printf{”Input  the  guassian  window  xmean  and  xvariance 
scanf(”%f  %{”,  &xmean,  &xvar); 
printf(”%f  %{\n”,  xmean,  xvar); 

printf(”Input  the  guassian  window  ymean  and  yvariance  :”); 
scan{(”%f  %{”,  &ymean,  &yvar); 
printf(”%f  %{\n”,  ymean,  yvar); 

♦♦♦♦</ 

/♦♦  Allocate  Memory  for  Arrays  **/ 
printf("Allocating  Memory  \n"); 
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outf ile\n\n"); 


p  =  matrix(l,  row,  1,  col); 
w  =  matrix(l,  row,  1,  col); 


printf("  raading  inage  of  Xd  row  and  Xd  col\n",  row,  col); 


for(j=l;  j<row;  j++) 
for(i=l;  i<col;  i++) 
f8canf(fin,  ”Xf\n'*,&p(i](j]); 


printf("  Calcualting  vlndov  \n“); 

normal  =  1/(2  ♦  pi  *  xvar  *  yvar); 

for(i=l;  i<row;  i++) 
for(j=l;  j<col;  j++){ 

w[i][j]  =  exp  ((double)  -0.5  ♦(  SQ((i— xmean)/xvar)+  SQ((j-ymean)/yvar))); 
p[i]U]  =  P[i]6]  *  w[i][j]; 

/**prmtf(”%d,  %d,  %g\a”,  i,j,  w[i][j]);**/ 

} 

re8cale(p,  row,  col); 
re8cale(w,  row,  col); 

{oi(j=l;  j<row;  j++) 
for(i=l;  i<col;  i++) 

fprmtf(fout,  "X4.0f\n  ",  p(i][j]); 

fo  =  fopen("Bind.dat",  "w"); 

for(j=l;  j<row;  j++) 
for(i=l;  i<col;  i++) 

fprintf(fo,  "X4.0f\n  ",  w[i][j]); 

fclo8e(fo); 

fclo8e{fout); 

fclose(fin); 

} 

void  rescale( output,  row,  col) 
float  -^^-output; 
int  row,  col; 

{ 

int  NEWJdAX,  NEW  JdIN,  i,  j,  count; 

float  min,  max; 

NEW_MAX  =  255; 

NEW_MIN  =0; 


/***  priatfC\n!!!  FINDING  MAX  **/ 
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!**  Check  for  the  max  and  min  value  in  the  data  **/ 
iiun=max=output[l][l]; 

coant=0; 


for(j=l;  j<row;  j++) 

for(i=l;  i<col;  i++){ 

if(output  [i]  [j]  >  max){ 

/♦*  pnat{(”%^n”,  max);  **/ 
max=output  [i]  [j] ; 

if(oatput[l]p]<min){ 

min=output[i][j]; 

/**  pnatf(”%l\a”,  min);**! 

} 

count++; 

} 

/♦♦*  print/]f”\n  SAMPLES  =  %d\n”,  count); 
printf(”  max  =  %{ min  =  %f\n”,  max,  min);  ♦♦V 


/***  Now  translate  data  and  write  to  output  file  ***/ 
for(j=l;  j<row;  j++) 

for(i=l;  i<col;  i++){ 

output[i][j]  =  ((output[i][j]— min)*(NEW_MAX— NEW^IN)/(max— min)  +  NEW3IIN); 

} 

/****  printf(”!!!  All  Done  !!!\n\n”);  ♦♦♦*/ 

} 


CJateS.c 

NAME:  CJateS.c 

INVOKED:  CJateSm  image!  image2  outiile 
DATE:  July  1991 

DESCRIPTION:  This  program  centers  image2  on  image!  and  places  the  centered 
iinage2  in  outiile. 

SUBROUTINES  CALLED:  CorrelateQ,  maxJindQ,  CJate3(),  shiftQ,  rescaleQ 
FUTURE  MODIFICATIONS/BUGS:  none 

t********************************************************************/ 

^include  <8tdio.h> 

^include  <mat.h.h> 

#define  skipJine  gets(junk) 

#define  IDX(a,b,c)  (a)*(c)+(b) 
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#define  SQR(a)  (a)*(a) 

#deime  loopi(A)  for(i=0d<(A);i++) 

#define  loopj(A)  for(j=0u<(A);j++) 

#defuie  loopij(A,B)  for  (i=0;  i<(A);  i++)\ 
for  (j=0;  j<(B);  j++); 

float  ♦vector(); 
void  foum(); 

/***  void  dodipO;  ***{ 
void  Correlate(); 
void  maxJind(); 

void  CJate3(); 

void  shift(); 

void  re8cale(); 

main(argc,argv) 
int  argc; 
char  *argv[|; 

{ 

int  x,y4  J,  Row,  Col, sub-row,  sub-col,  downsample; 

FILE  *fout,  ^datJilel,  «datJile2,  *outJile; 

char  filename[81],  in-string[81],  junk[256]; 

float  *xl,  ♦x2,  *vector(),  ♦♦matrix(),  ♦♦ims^el,  ««image2; 

signed  int  location[3]; 


if  (argc  /  4)  { 

printf("!!!  The  conmand  line  should  be  !!!:\n\n  C.lateS.n  image!  image2  outfile 
\n\n"); 
exit(O); 

} 

l^^j^^^t************  Set  Up  Files  ♦♦♦***♦♦****♦*♦♦♦*♦♦♦♦*</ 

if  ((dat-fllel  =  fopen(argv[l],  "r"))  ==  NULL)  { 
printf("I  can't  open  the  image!  file"); 
exit(— 1); 

} 

if  ((dat.flle2  =  fopen(argv[2],  "r"))  ==  NULL)  { 
printf("I  can't  open  the  output  file"); 
exit(-l); 

} 

if  ((font  =  fopen(argv[3],  "w"))  ==  NULL)  { 
printf("I  can't  open  the  image2  file”); 
exit(-l); 

} 

l*****pTiat{(”Input  the  number  of  rows  and  cols  in  Row  Col 
scanf(”%d  %d”,  &Row,  &CoI); 
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Row  =  128; 

Col  =  128; 

I*  dynamically  allocate  memory  */ 

imagel  =  matrix(0,  Row-1,  0,  Col-1); 

image2  =  matrix(0,  Row-1,  0,  Col-1); 

/*  Initialize  matrix  **/ 

for(y=0;  y<Row;  y++) 
for(x=0;  x<Col;  x++){ 

imagsl[x][y]=  image2[x][y]  =  0.0; 

} 

read  the  training  faces  in  as  columns  ****/ 
prmtf("M!  Opening  and  Reading  Ref  Face  l:\n  "); 

for(y=0;  y<Row;  y++) 
for(x=0;  x<Col;  x++){ 

f8canf(datJllel,  "Xf\n",&imagel[x][y]); 

} 

fclo8e(dat_filel); 

printfC'!!!  Opening  and  Reading  Face  to  be  Centered: \n  "); 
for(y=0;  y<Row;  y++) 
for(x=0;  x<Col;  x+4-){ 

f8canf(dat_file2,  ''Xf\n'',&image2[x][yJ); 

} 

fclo8e(dat_file2); 


START  FACE  CENTERING  *♦**/ 

downsample  =  1;  /***  downsample  allows  for  faster  centering  by  desampling  the  image 

sub  jow  =  Row/downsample;  /****  down  sample  image  for  faster  correlation  **f 
sub. col  =  Col/downsample; 

printf("init  vectoreNn"); 

xl  =  vector(0,2*subjow+sub_col-l); 
x2  =  vector(0,2*8ubj'ow*sub_col-l); 


prmtf(" initialize  vectorsXn"); 


loopi(8ubj‘ow*8ubjcol){ 
xlf2*i]  =  xl[2*i+l]  = 
2*il  =  x2[2*i+lj  = 


xl 

x2 


} 


0.0; 

0.0; 


/*t********************  PUT  MATRIX  INTO  VECTOR  FOR  FFT  OPS  *******/ 
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loopi(8ub_row)  { 

loopj(8ub-col)  (float)  xl[2*(i*8ubxol+j)]  =  imagel[(j+downsample)][(i*downsample)]; 
loopj(8ub^ol)  (float)  x2[2«(i*subjcol+j)]  =  image2[j*downsample][i4>downsample]; 

} 

CJa,te3  return  the  amount  of  borizontal(i)  or  verticalQ)  shift  needed  ****f 
C  Jate3(xl,  x2,  8ub  jow,  6ub_col,  location); 

prmtf("  i  location  Xd  j  location  Xd  \n",  location[0],  locationfl]); 

/***Some  Error  Correction  ****j 

if(location(0]>(8ub_col/2))  locationfO]  =  — (8ubjcol  —  location[0]); 
if(location[l  >(8ubjrow/2))  locationfl]  =  —(sub-row  —  locationflj); 
if(location[0  <(sub_col/2))  location[0]  =  locationfO]; 
if(location[l]<(8ub.row/2))  locationfl]  =  locationfl]; 

location[0]  =  location[0]  *  downsample; 
locationfl]  =  location[l]  *  downsample; 


4<«««4<««i»*4<********4<««*4^*<*«  Actually  Shift  Image  with  shift  ****/ 

shift(image2,  Row,  Col,  location); 

fi,iL**it*<^^it**********^'*********  Make  0  to  255 
rescale(image2.  Row,  Col); 

for  (y  =  0;  y  <  Row;  y++) 
for  (x  =  0;  X  <  Col;  x++) 

fprintf(fout,  "X4.0f\n  ",image2[x](y]); 
fclose(fout);  ♦*♦*/ 

printf("  END  OF  MAIN\n"); 

} 

f^^t*************************************************************** 
Subroutine  CJate3 

INVOKED:  CJate3(xl,  x2,  Row,  Col,  location) 

DATE:  July  1991 

FUTURE  MODIFICATIONS/BUGS:  none 

^^f*:^^:t,ti^*^t**********************************************************/ 

void  CJate3(xl,  x2.  Row,  Col,  location) 
int  Row,  Col; 
signed  int  location  Q; 
float  xlQ,  x2Q; 

{ 

FILE  ♦outJile; 
int  n(2],  i,  j; 
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float  ^output,  temp,  <»*matrix(),  *tmatJoatput,  *vector(); 


/♦*  Allocate  Memory  for  Arrays  **/ 

/***  pr}nt{(” Allocating  Memory  \n”);  ****/ 

output  =  vector(0,2*Row*Col-l); 
matjoutput  =  matrix(0,  Row— 1,  0,  Col— 1); 

/**  Assign  Initial  Array  Values  **/ 

n[0]  =  Col; 
n[lj  =  Row; 


printf(  "Corellat  ing\n"  ); 
Correlate(xl,  x2,  output,  n.  Row,  Col); 


printf( "Storing  resultsXn"); 

/**  Store  The  Magnitude  Results  ♦*/ 
loopi(Row) 
loopj(Col)  { 

temp=sqrt((double)SQR(output[2*(i*Col+j)])  +(double)SQR(output[2*(i*Col+j)+l])); 
/**temp=temp/sqrt((double)SQR(outputfOj)  i-(doub!e)SQR(output[l]));**/ 

/***  fprmtf(out_fiie,”%4.2f\n”,  temp); 
matjDutput[j][i]  =  (float)  temp; 

} 

maxJind(matjoutput,  Row,  Col,  location); 

/♦♦♦  fclose(outJile);  ****/ 

}  A*  *V 


Subroutine  Correlate 

INVOKED:  Correlate(inputl,  input2,  output,  n,  Row,  Col) 

DATE:  July  1991 

FUTURE  MODIFICATIONS/BUGS:  none 

^^t1,^^,*^,:^^^^t^t^t*******************************************************/ 

void  Correlate(inputl,  input2,  output,  n.  Row,  Col) 
float  inputl[],input2Q,  outputQ; 
int  nQ,  Row,  Col; 

{ 

int  i; 

float  *templ,  «temp2; 

tempi  =  vector(0,2*Row*Col-l); 
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temp2  =  vector(0,2*Row*Col— 1); 


loopi(2*Row*Col){ 
templ[i]  =  inputl[i]; 
temp2[i]  =  input2[ij; 

} 

/**  TBike  JFburier  Transform  of  Input  Functions  **/ 

founi(templ— 1,  n— 1,  2,  1); 
fourii(temp2— 1,  n— 1,  2,  1); 

/**  Conjugate  One  of  The  Fourier  Transforms  **/ 
loopi(Rx>w*Col) 

temp2[2*i+i]  =  — temp2[2*i+l]; 

/**  Multiply  Fourier  Transforms  Together  **/ 

loopi(Row*Col){ 

/**  Real  Component  **/ 

output[2«i]  =  templ[2*i]«temp2[2«i]  —  templ[2«+l]*temp2[2*i+l]; 

!**  Imaginary  Component  **/ 

output[2*i+l]  =  templ[2*i]*temp2[2*i+l]  +  temp2[2*i]*templ[2*i+l]; 

} 

/**  Take  Inverse  Transform  to  obtain  Correlation  **/ 

fourn(output— 1,  n— 1,  2,  —1); 

/*♦  Rescale  to  get  proper  magnitude  **/ 

loopi(2*Ilow*Col) 
output[i]  /=  Row*Col; 

/♦♦*♦♦♦**♦♦*****♦*♦♦♦♦♦♦♦♦♦*****♦♦**♦*♦♦*♦**♦******♦**********♦****♦**♦*+** 
The  result  of  the  correlation  is  that  first  element  of  the  output 
matrix  is  for  zero  shift,  the  next  element  for  shift  one  to  the  right  and 
so  on.  This  puts  results  into  a  format  which  humans  can  understand. 
****0*********************************************************************^ 


!**  Free  up  the  memory  when  finished  ♦♦/ 

free-vector(templ  ,0,2*Row*Col- 1); 
free.vector(temp2,0,2*Row*Col- 1); 

}  /♦♦♦♦♦♦♦♦♦♦♦*♦***♦♦♦♦  End  of  Correlate  Routine  *********************^/ 
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Subroutine  maxJiDd  finds  maximum  value  in  a  matrix  and  returns  location 
INVOKED:  maxJind(matjoutput,  row,  col,  location) 

DATE:  July  1991 

FUTURE  MODIFICATIONS/BUGS:  none 

^:^:****^^^i********^t**************************************************^ 


void  maxJind(matjoutput,  row,  col,  location) 
float  **matjoutput; 
int  row,  col; 
signed  int  locationQ; 

{ 

int  i,  j,  count,  tempJ,  temp^; 
float  max; 

/***printf(”\n!!!  FINDING  MAX  ///\n\n”>*V 

/**  Check  for  the  max  and  min  value  in  the  data  **/ 

max=matjoutput[0][0]; 

count=0; 

for(j=0;  j<row;  j++) 
for(i=0;  i<col;  i++){ 

if(matjoutput[i][j]>max)  { 

max=matjoutput[i][j]; 
tempJ=i;  temp_j=j; 

} 


count++; 

} 

location[0]=tempJ;  location[l]=temp_j; 
if(location[0] >col)  location[0]=0; 
if(location[lj>row)  location[l]=0; 

/*♦♦♦ 

printfr\n\n  SAMPLES  =  %d\n”,  count); 
printf(”  max  =  %f\n”,  max); 

printf(”i  location  =  %d  j  location  =  %d  \n”,  locationfO],  locationfl]); 
printf(”\n!!!  All  Done  !!!\n\n”); 

} 

void  shift(image.  Row,  Col,  location) 
float  *4>image; 
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int  Row,  Col; 
signed  int  locationQ; 

{ 

/*t***t**^f^L*^^****t******************************************************** 

NAME:  shift 

DESCRIPTION:  This  routine  takes  a  image  shifts 

*****************************************t****t*************************^ 

int  f^tempJmage,  i,x,y,  xshift,  yshift,  abs(); 
float  **matrix(),  ^^sliiftedJmage; 
int  n,k,  new  jow,  new_col; 

xshift  =  location[0]; 
yshift  =  location[l]; 

printf("  xshift  ■  Xd  yshift  ■  Xd  \n'',  xshift,  yshift); 

printf("  Allocating  Memory \n"); 

newjow=  (int)  Row+  2  ♦  abs(yshift); 
new_col=(int)  Col+  2  *  abs(xshift); 

printf("Xd  Xd  Xd  Xd\n",  -abs(xshift),  new.col,  — abs(yshift),  newjow); 
shiftedJmage  =  matrix(  — abs(xshift),  new_col,  — abs(yshift),  newjow); 

printf("  initialize  matrixXn"); 

/***  initialize  matrix  ***/ 

for(y=  -  abs(yshift);  y<newjow;  y++) 
for(x=  —  ab8(x8Uft);  x<  new_col;  x++) 
shifted  Jmage[x][y]= 127.0; 

printf("  shifting  imageNn"); 

for(y=0;  y<Row;  y++) 
for(x=0;  x<Col;  x++){ 

8hiftedJmage[x+xshift][y+yshift]  =  image[x][y]; 


for  (y  =  0;  y  <  Row;  y++) 
for  (x  =  0;  X  <  Col;  x++) 

image[x]  [y] =shifted  Jmage[x]  [y] ; 


} 


NAME:  rescale 

DESCRIPTION:  This  routine  takes  a  rescale  an  range  to  0  to  255 
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void  rescale(output,  row,  col) 
float  **output; 
int  row,  col; 

{ 

int  NEW_MAX,  NEW J^IN,  i,  j,  count; 

float  min,  max; 

NEW-MAX  =  255; 

NEW-MIN  =0; 


printf("\n! ! !  FINDING  MAX  !!!\n\n''); 

/**  Check  for  the  max  and  min  value  in  the  data  **/ 
min=max=output[0][0]; 

count =0; 


for(j=0;  j<row;  j++) 

for(i=0;  i<col;  i++){ 

if(out  put  [i]  [j]  >  max)  { 

/*  prmtf(”%^a”,  max);  */ 
max=output  [i]  [j] ; 

} 

if(output  [i]  p]  <  min)  { 

min=output[i][j]; 

/**  prmtf(”%l\a”,  min);**! 
count++; 


/**  priat{(”\n  SAMPLES  =  %d\n”,  count); 
pnntf(”  max  =  %{ min  =  %f\n”,  max,  min);  ♦♦♦*/ 


/***  Now  translate  data  and  write  to  output  hie  ***/ 
for(j=0;  j<row;  j4-+) 

for(i=0;  i<col;  i++){ 

output[i][j]  =  ((output[i][j]-min)*(NEW.MAX— NEW_MlN)/(max-min)  +  NEW_MIN); 


} 


fft-trunc.c 

/****iL^*t^i^tt*****************************^*********************** 

NAME:  fftJrunc.c 

INVOKED:  fft.trunc  inhle  outhle  imsize  order 
DATE:  June  1991 
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DESCRIPTION:  This  program  takes  the  FFT  of  an  image  and  if  order  is  greater  than 
zero  it  provides  the  magnitude  of  the  FFT  upto  order. 

SUBROUTINES  CALLED:  truncate,  fourn,  doSip 
FUTURE  MODIFICATIONS/BUGS:  Ugly  code 

4,^,^t^i^I^^^^:^^,^,^l^l^c^L:^:^t************************************************>^ 

#include  <device.h> 

^include  <stdio.h> 

^include  <math.h> 

#define  SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr 

#define  loopi(A)  for(i=0;i<(A);i++) 

#define  loopj(A)  for(j=:0y<(A);j++) 

^define  loopij(A,B)  for  (i=0;  i<(A);  i++)\ 
for  (j=0;  j<(B);  j++); 

#define  SQ(A)  (A*A) 

char  OQtname[80]; 
subroutine:  truncate 

^^t***************************************************^f*>(f^ 

void  tnincate( data, size, order, out jdata, name) 
float  data[],outjdata[]; 
int  8ize,order; 
char  nameQ; 

/**  truncate  an  FFT  to  order  desired  ***/ 

FILE  *fout; 

int  ij,  p,  rowjoffset,  coljoffset,  Row,Col; 
float  tempr; 

Row=8ize; 

Col=size; 

p=0; 

rowjoffset  =  (int)  Row/2-1; 
coljoflfset  =  (int)  Col/^ 

/♦  pritttf(”%s\n”,  name);  */ 
printf("Xf\n",data[8128]); 
fout=fopen(outname,  "a"); 

for(j=(-l*order);j<Oy++) 

{ 

for(i=(-l+order);i<order;i++) 

{ 

out^ata[p]  =  data[(colj}ffset— i)+((rowj3ffset— j)+Col)]; 
fprintf(fout,"X. loi  ",out.data[p]); 

P++; 
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} 


} 

fprmtf(fout,  "XsXn",  name); 

fprintf(fout , "  \n"  ) ; 
fdo8e(fout); 

} 

subroutine:doSip  -  Sips  the  fourn  fFt 

******************************************llL^ft***********^ 

void  dofiip(data,8ize) 
float  dataQ; 
int  size; 

{ 

/***convert  fount  format  to  format  normal  humans  can  use  **/ 
int  i  j,  rowjoflPset,  coljoffset,  Row, Col; 
float  tempr; 

Row=8ize; 

Col=size; 

rowjoflPset  =  (int)  Row/2  ♦  Col; 
coljoffset  =  (int)  Col/2; 

for( j=0u  <  (Row*  Col/2)  y + = Col) 

for(i=0;i<Row;i++)  /*  top  half  to  botttom  half  swazp  */ 

tempr  =  data[i+j]; 

data[i+j]  =  data[i+j+rowjofFset]; 

data(i+j+rowjoff8et]  =  tempr; 

} 

for(j=0;  j<((Row-l)*Col);j+=Col) 
for(i=0;  i<Col;  i++)  /*  left  half  to  right  half  swap  */ 

{ 

tempr  =  data[i+j]; 

data[i+j]  =  data[i+j+coLoffset]; 

data[i+j+coljjffset]  =  tempr; 

} 


} 

subroutine: fount  -  Numerical  Rec  in  C  FFT  routine 

4iiti%«*«^«4i«t*********************************************y 


void  foum(data,nn,ndim4sign) 
float  dataQ; 
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int  nn[],ndim4sign; 

{ 

int  il424342rev43rev4pl4p24p3?ifpl4fp2; 
int  ibit4dim4cl,k2^,nprev,nrem,ntot; 
float  tempi, tempr; 
double  theta, wi,wpi,wpr,wr,wtemp; 

ntot=l; 

for  (idim=l;i(iim<ndim;idim++) 
ntot  *=  nn[idim]; 
nprev=l; 

for  (idim=ndim;idim>14dim — )  { 
n=nn[idim]; 
nrem=:ntot  /  (n^nprev); 
ipl=nprev  <  1; 
ip2=ipl*n; 
ip3=ip2*nrem; 
i2r6v=l* 

for  (i2=l;i2<ip242+=ipl)  { 
if  (i2  <  i2rev)  { 

for  (il=i2;il<i2+ipl-2;il+=2)  { 
for  (i3=il43<ip3,i3+=ip2)  { 
i3rev=i2rev+i3-i2; 
SWAP(data[i3],data[i3rev]); 
SWAP(data[i3+l],data[i3rev+l]); 


} 

ibit=:ip2  >  1; 

while  (ibit  >  ipl  &&  i2rev  >  ibit)  { 
i2rev  — =  ibit; 
ibit  >=  1; 

} 

i2rev  +=  ibit; 

} 

ifpl=ipl; 

while  (ifpl  <  ip2)  { 
ifp2=ifpl  <  1; 

theta=isip+6.28318530717959/Cifp2/ipl); 

wtemp=8in(0.5+theta); 

wpr  =  — 2.0*wtemp*wtemp; 

wpi=8in(theta); 

wr=1.0; 

wi=0.0; 

for  (i3=l;i3<ifpl;i3+=ipl)  { 

for  (il=i3;il<i3+ipl-2;il+=2)  { 
for  (i2=il;i2<ip3;i2+=ifp2)  { 
kl=i2; 
k2=kl+ifpl; 

tempr= wr*data[k2]  -  wi+data{k2+ 1]; 
tempi=wr*data[k2+l]+wi4>data{k2]; 
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} 


} 


k2]=data[kl]— tempr; 


data[ 

data[k2+ 1] =data[k  1 + 1] -tempi; 
dataf 
data 


kl]  +=  tempr; 
kl+1]  +=  tempi; 


wr=(wtemp=wr)*wpr— wi*wpi+wr; 
wi=wi*wpr+wtemp*wpi+wi; 

} 

ifpl=ifp2; 

} 

nprev  *=  n; 


/*♦********♦♦♦♦♦**♦**♦♦*♦**♦****♦*♦♦♦♦**♦♦♦♦*★*♦♦♦♦*♦♦♦♦♦♦ 

MAIN 


^f^f^^t^t****************************^*********************’^ 

main(argc,argv) 
int  argc; 
char  *argv[]; 

{  FILE  ♦fin,  ♦font; 

float  ♦output, ♦input, ♦truncjout; 
float  norm; 
float  ♦vector(); 
void  ♦free.vectorO; 
char  name[30]; 

int  ij,  nn[i],  ndim,  isign,  new  .order,  order,  image.size; 
if(argc  ^  5){ 


printf("!!!  The  command  line  should  be  !!!:\n\n  fft.tnmc  infile  outfile  imsize 
order  \n\n"); 
exit(O); 

} 


imagejsize=atoi(argv[3]); 

order=atoi(argv[4]); 

sprintf(outname,  "XsXd","netf ft .ng",order); 

/^■♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦set  up  dynamic  allocation****************^ 
input  =  vector(0,2^imagejsize^imagejsize— 1); 
output  =  vector(04niage.size^image.8ize-l); 


/*^i*t**************  Set  Up  Files  ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦^ 
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if  ((fin=fopen(axgv[l],"r"))  ==  NULL){ 
prmtf("I  can't  open  the  input  file"); 
exit(— 1); 

} 


if  ((fout=fopen(argv[2],"w"))  ==  NULL){ 
printf("I  can't  open  the  output  file"); 
exit(-l); 

} 

/**************Rje&d  File  *****0*************^ 

loopi(2*imagejBize<t>image^ize— 1)  /*  initialize  array  to  zero  */ 
input[i]  =  0.0; 

loopi(imagej5ize*image^ize— 1)  l*read  dataia  the  fourn  format 
{  j*  see  numerical  recipes  in  c  ♦/ 

f8canf(fin,  "Xf  \n",  &input[i*2]); 

} 


fclo8€(fin);  f*cIose  input  file  */ 


/*  Initialization  parameters  for  FFT  */ 

nn[0]=:iinagejsize;  /♦  size  of  mput  lAW  fournQ  */ 
nn[l]=imagejBize; 

ndim=2;  /*  two  dim  FFT  ♦/ 

i8ign=l;  /*  FFT  */ 

fourn(input— 1,  nn— 1,  2,  1); 


/*♦**♦♦♦♦♦♦*  Find  Fourier  Magnitude  ♦♦♦♦******y 
j=0; 

for(i=0;i<(2«image.size«image-8ize-l);  i+=2) 

{ 


output[j]=8qrt(  (double)  SQ(input[i])  + 

(double)  SQ(input[i+l])); 

j++; 

} 

norm=output[0];  /*  d.c  component  used  for  normalization  ***/ 
printf("X4 .  Of  \n", norm); 
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/*****  doAip*******************************/ 


doflip(oatput4niage^ize);  j*  converts  foura  format  to  human  format  */ 
printf("X4 .41\n"  ,output[8128]); 

loopi(image  jBize*image  jsize)  { 
output[i]=  1000*oatput[i]/norm; 
fprintf(fout,  "Xf\n",  output(i]); 

} 

/*****  normalize  and  write  output  of  FFT  in  argv[2]  hie  *****/ 

^t*****  truncate  tt^*********************************************************/ 
I*  truncate  takes  fft(output)  of  size(image^ize)  and  truncates  the  FFT  to  **/ 
f*  order  specified  plus  d.c.  the  array  is  returned  in  trunc.out,  the  argv[2]*/ 

/*  is  used  as  a  heaider  when  truncate  writes  the  output  in  netfft.dat 


if(or<ler  ^  0){ 

newjorder  =  2*order+l; 

truncjout  =  vector(04magejsize*image-size— 1); 

t ruiicate(output  4mage^ize, order ,trunc^ut ,  argv  [2] ) ; 


/♦  free.vect  or(t  run  Cjou  t  ,0,image_size*  image-size-  J  );*/ 

) 


free.vector(input,0,2#image^ize*image-size-l); 
free_vector(output  ,04inage^ize«image_size— 1 ); 
fclose(fout); 


} 


noise.c 

#mclude  <8tdio.h> 
#include  <math.h> 
#define  Ml  259200 
#define  lAl  7141 
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#define  ICl  54773 
#defuie  RMl  (1.0/Ml) 

#deiine  M2  134456 
#define  IA2  8121 
#define  IC2  28411 
#define  RM2  (1.0/M2) 

#define  M3  243000 
#deiiiie  IA3  4561 
#deiine  ICS  51349 

main(argc,argv) 
int  argc; 
char  ♦argvQ; 

FILE  *datJllel,  #datJiie2; 
float  gasdevO; 
float  x4nval,  outval; 
int  i,var,  dum; 

if  (argc  #  5)  { 

printf("!M  The  command  line  should  be  !!!:\n\n  noise  inimagel  outimage  seed 
var\n\n''); 
exit(O); 

} 

printf("Input  the  reandom  seed  and  var 
8canf("Xd  %d",  &dum,  &vaj); 
dum  =  atoi(arsv[3]); 
var  =  atoi(argv[4]); 

if  ((datJilel  =  fopen(argv[l],  "r"))  ==  NULL) 

{ 

printf("I  can't  open  the  imagel  file"); 
exit(-l); 

} 

if  ((dat-flle2  =  fopen(argv[2],  "w"))  ==  NULL) 

printf("I  can't  open  the  output  file"); 
exit(-l); 

} 


while(fscanf(dat_fllel,  "Xf\n",  &inval)  ^  EOF){ 


outval  =  inval  +  (var  ♦  (float)  gasdev(  (int)&dum)); 
fprmtf(datJiIe2,  "X4.0f\n",  outval); 

} 


} 
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float  gasdev(idum) 
int  *idum; 

{ 

static  int  iset=:0; 
static  float  gset; 
float  fac,r,vl,v2; 
float  ranlQ; 

if  (iset  ==  0)  { 
do  { 

vl=2.0*ranl(idum)— 1.0; 
v2=2.0*ranl(idum)— 1.0; 
r=vl*vl+v2*v2; 

}  while  (r  >  1.0); 

fac=sqrt(-2.0*log(r)/r); 

gset=vl*fac; 

iset=l; 

return  v2*fac; 

}  else  { 

iset=0; 
return  gset; 

} 

} 

float  ranl(idam) 
int  «idum; 

{ 

static  long  ixl,ix24x3; 
static  float  r[98]; 
float  temp; 
static  int  ifl'=0; 
int  j; 

void  nrerror(); 

if  (*idum  <  0  II  iff  ==  0)  { 
iff=l; 

ixl=(ICl-(*idum))  %  Ml; 
ixl=(IAl*ixl+ICl)  %  Ml; 
ix2=ixl  %  M2; 
ixl=(IAl*ixl+ICl)  %  Ml; 
ix3=ixl  %  M3; 
for  (j=lu<97u++)  { 

ixl=(IAl*ixl+ICl)  %  Ml; 
ix2=(IA2*ix2+IC2)  %  M2; 
r[j]=(ixl+ix2*RM2)*RMl; 

} 

«idum=l; 

} 

ixl=(IAl*ixl+ICl)  %  Ml; 
ix2=(IA2*ix2+IC2)  %  M2; 
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ix3=(IA3*ix3+IC3)  %  M3; 
j=l  +  ((97*ix3)/M3); 

if  (j  >  97  II  j  <  1)  nrerror("RANl:  This  cannot  happen."); 
temp=rp]; 

r(j]=(ixl+ix2*RM2)*RMl; 
return  temp; 

} 

#andef  Ml 
#undef  lAl 
#undef  ICl 
^undef  RMl 
^undef  M2 
^undef  IA2 
#undef  IC2 
#undef  RM2 
#undef  M3 
#undef  IA3 
#undef  IC3 


angle.c 

/^t*************************************************************** 

NAME:  angle.c 

INVOKED:  angle  weigbtoutfile  totalsize  numberofvectors 
DATE:  24  June  1991 

DESCRIPTION:  This  program  calculates  the  angles  between  eigenfaces  1  ~  8 
SUBROUTINES  CALLED:  none 

FUTURE  MODIFICATIONS/BUGS:  change  file  i/o  method  to  be  compatible  with 
recon.c 

t*******************************************************************’^ 

^include  <stdio.h> 

#mclude  <math.h> 

^define  SQ(A)  (A+A) 
main(argc,argv) 
int  argc; 
char  ♦argvQ; 

{ 

FILE  ♦facel,  4>face2,  «face3,  ♦face4,  ♦faceS,  ♦faced,  ♦face7,  ♦faceS,  ♦fout; 

FILE  ♦face^vg,  ♦fevex,  ♦feval; 
int  ij,  N,  k,  M,  nrot; 

float  ♦♦matrix(),  ♦vector(),  ♦♦A,  **Wj  ♦m^; 
float  ♦d,  temp; 

void  free_vector(),  free-matrix(),  eigsrt(),  jacobi(); 
double  sqrt(),  aco8(); 
char  atoi(); 

if(argc  ^  4){ 
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printf('’!!!  The  command  line  should  be  !!!;\n\n  angle  eeightoutf ile  totalsize 
numberofvectors\n\n"); 
exit(O); 

} 

/^:**^i*^^,ift*********  Set  Up  Files  ♦**♦♦*♦♦♦*♦♦♦*♦*♦♦***♦♦■»/ 
if  ((fout=fopen(argv[l],"w"))  ==  NULL){ 
prmtf("I  can't  open  the  angle  output  file"); 
exit(-l); 

} 

N  =  atoi(argv[2]); 

M  =  atoi(argv[3]); 


/*  dynamicaJiy  allocate  memory  */ 

A  =  matrix(l,r,l,M); 
w  =  matrix  (1,M,  1,M); 
mag  =  vector(l,  M); 

/•>♦♦♦  initalize  matrix  and  vectors  ****/ 

for(j=l;j<M;j++) 

for(i=l;i<N;i++) 

A[i][j]=0.0; 


for(k=l;  k<M;  k++) 
for(j=l;  j<M;  j+-r) 
w[j][k]  =  mag(k]  =  0.0; 


fiHi**^*  read  tbe  eigen  faces  in  as  columns  ****j 

printf("!!!  Opening  and  Reading  EigenFace  l:\n  "); 
facel=fopen("eigenf aceOl  .dat",  "r"); 
for(j=lu<Nu++){ 
f8canf(facel,"Xf\n",&A[j][l]); 

/*printf(”%4.0{\n”,  A[j][l]);^ 

} 

fclose(facel); 


printf("!!!  Opening  and  Reading  Eigenface  2:\n  "); 
face2=fopen("eigenf  ace02 . dat",  "r"); 
for(j=lu<Ny++){ 

f8canf(face2,  "Xf\n",  &A[j][2]); 

} 

fclose(face2); 


printf("!!!  Opening  and  Reading  Eigenface  3:\n"); 
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face3=fopen(''eigenface03.dat",''r"); 

for(j=lu<Ny++){ 

fscanf(face3,  "Xf\n",  &A(j][3]); 

} 

fclose(face3); 


printf("!!!  Opening  and  Reading  Eigenface  4:\n  "); 
face4=fopen("aigenf ace04.dat",  “r"); 

for(j=l-d<Nu++){ 

fscanf(face4,  "Xf\n",  &A[j][4]); 


} 

fclose(face4); 


printf("!!!  Opening  and  Reading  Eigenface  5:\n"); 

face5=fopen("eigenf aceOS . dat",  “r"); 
for(j=lu<Nu++){ 

fscanf(face5,  "Xf\n",  &A[j][5]); 


} 

fclo8e(face5); 


printf("!!!  Opening  and  Reading  Eigenface  6:\n  "); 

{'ace6=:fopen("eigenf ace06 .dat",  "r"); 

for(j=lij<Nu++){ 

f8canf(face6,  "Xf\n",  &A[j][6]); 


} 

fclo8e(face6); 


printf("!!!  Opening  and  Reading  Eigenface  7:  \n"); 

face7=fopen("eigenf ace07 .dat",  "r"); 

for(j=lU<Nu++)  { 

f8canf(face7,  "Xf\n",  &A[j][7]); 

} 

fclo8e(face7); 


printf("!!!  Opening  and  Reading  Eigenface  8:\n  "); 

face8=fopen("eigenf  aceOS .  dat" ,  "r"); 
for(j=lU^Na++)  { 

fscanf(face8,  "Xf\n",  &A[j][8]); 


} 

fclose(face8); 
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printf(“!!!  Noxnializlng  :\n  "); 

/***  Normalizing  Data  by  dividing  by  255  ***/ 


for(j=ly<M;j++) 

for(i=l;i<N;i++) 

A[i]D]=A[i]Dl/255; 


printf("!!!  Calculating  magnitudesW,  M); 
for(j=l;  j<M;  j++){ 
for(i=l;  KN;  i++) 

=  A[i][j]*  A[i][j]+  mag[j]; 
mag[j]  =  8qrt(  (double)  mag[j]  ); 

/***  print{(”  Wmt!!\n”);  ***</ 

} 

printf("!!!  Calculate  dot  products  \n",  M); 
for(k=l;  k<M;  k++){ 
for(j=l;  j<M;  j++) 
for(i=l;  i<N;  i++) 

w[j][k]  =  A(i][j]  *  A[i][k]  +  wLj][k]; 

} 

printf("!!!  Cadculate  cos  alpha  \n",  M); 
for(k=l;  k<M;  k++) 
for(j=l;  j<M;  j++) 

w[j][k]  =  wfi][k]/(mag[j]*mag[k]); 

printf("!!!  Write  cos  alpha  \n",  M); 

for(k=l;  k<M;  k++) 
for(j=l;  j<M;  j++) 

fprintf(fout, "Angle  in  degrees  between  uXd  and  uXd  is:  X3.2f\n", j, k,57.2957*acos( 
(double)  w[j][k])); 

fclose(fout); 

printfC'!!!  ALL  DONE  \n",  M); 

} 


kl-med 

f*****^i***>li********:^**************1f*******************iHt********** 

NAME:  fftJrunc.c 

INVOKED:  flt.trunc  inHIe  outlile  imsize  order 
DATE:  June  1991 

DESCRIPTION:  This  program  takes  the  FFT  of  an  image  and  if  order  is  greater  than 
zero  it  provides  the  magnitude  of  the  FFT  upto  order. 

SUBROUTINE'S  CALLED:  truncate,  fourn,  doSip 
FUTURE  MODIFICATIONS/BUGS:  Ugly  code 
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#include  <device.h> 

^include  <stdio.h> 

^include  <math.h> 

#define  SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr 

#define  loopi(A)  for(i=0;i<(A);i++) 

#defiiie  loopj(A)  for(j=0;j<(A)y++) 

^define  loopij(A,B)  for  (i=0;  i<(A);  i+-l-)\ 
for  (j=0;  j<(B);  j++); 

#define  SQ(A)  (A*A) 

char  outname[80]; 

f^ifi^L^L**^*^i**^^f^n^iHr^^Ht^L************************************ 

subroutine:  truncate 

^t^^,^,^^^^:^*^t^L*^i:^^^^*^,:^^,^**^Ht^t^^^^^:^^^t************************/ 

void  tnincate(data, size, order, out jdata,name) 
float  data[],outjdata[]; 
int  size,order; 
char  nameQ; 

{ 

/**  truncate  an  FFT  to  order  desired  ***/ 

FILE  ♦fout; 

int  ij,  p,  rowjoflfset,  coLoffset,  Ilow,CoI; 
float  tempr; 

Row=8ize; 

Col=8ize; 

p=0; 

rowjoffset  =  (int)  Row/2-1; 
coljoffset  =  (int)  Co\/% 

/*  printf(”%s\n”,  name);  */ 
printf(  "Xf  \n"  ,data[8128]); 
fout=fopen(outnaine,  "a"); 

for(j=(-l+order)y<Oy++) 

{ 

for  (i= ( — 1  ♦order)  ;i  <order  ;i + + ) 

{ 

outjdata[p]  =  data[(coljoffset— i)+((rowjoffset— j)+Col)]; 
fprintf(fout,''X. iOf  ”,out^ata[p]); 

P++; 

} 

} 

fprintf(fout,  "XaXn",  name); 
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fprintf(fout,"\n"); 

fclo8e(fout); 

} 

subroutine:do3ip  -  Sips  the  fount  fft 

4,i,i,^:^i,i,4**t*********************************************^/ 

void  doflip(data,size) 
float  dataQ; 
int  size; 

/***coavert  fourn  format  to  format  normal  humans  can  use  **/ 
int  ij,  rowjoffset,  coljoflTset,  Row, Col; 
float  tempr; 

Row=size; 

Col=8ize; 

rowjoifset  =  (int)  Row/2  ♦  Col; 
coljofTset  =  (int)  CoXfi, 

for(j=0;j<(Row*Col/2);j+=Col) 

for(i=0;i<Row;i++)  /*  top  half  to  botttom  half  swazp  */ 

{ 

tempr  =  data[i+j]; 

data[i+j]  =  data[i+j+rowjofFset]; 

data[i+j+rowjoff8et]  =  tempr; 

} 

for(j=0;  j<((Row-l)*Col);j+=Col) 
for(i=0;  i<Col;  i++)  /*  left  half  to  right  half  swap  */ 

{ 

tempr  =  data[i+j]; 

data[i+j]  =  data[i+j+coljoffset]; 

data[i+j+coljofrset]  =  tempr; 


} 

subroutine:fourn  -  Numerical  Rec  in  C  FFT  routine 

t********************************************************/ 

void  fourn(data,nn,ndim48ign) 

float  dataQ; 

int  nnQ,ndim4sign; 

{ 

int  il424342rev43rev4pl4p24p34fpl4fp2; 
int  ibit  4dim  ,k  1  ,k2  ,n  ,npre v  ,nrem  ,ntot ; 
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float  tempi, tempr; 

double  theta, wi,wpi,wpr,wr,wtemp; 


ntot=l; 

for  (idim=l;idim<ndim;i<lim++) 
ntot  *=  nn[idim]; 
nprev=l; 

for  (idim=ndim;idim>l;idim — )  { 
n=nn[idim]; 
iirem=ntot  /  (n*nprev); 
ipl=nprev  <  1; 
ip2=ipl*n; 
ip3=ip24>nrem; 
i2rev=l; 

for  (i2=l;i2<ip2;i2+=ipl)  { 
if  (i2  <  i2rev)  { 

for  (il=i2;il<i2+ipl-2;il+=2)  { 
for  (i3=il;i3<ip3;i3+=ip2)  { 
i3rev=i2rev+i3-i2; 
SWAP(dataii3],data[i3rev]); 
SWAP(data[i3+l],data(i3rev+l]); 


} 

ibit=ip2  >  1; 

while  (ibit  >  ipl  &&  i2rev  >  ibit)  { 
i2rev  -=  ibit; 
ibit  >=  1; 

} 

i2rev  +=  ibit; 

} 

ifpl=ipl; 

while  (ifpl  <  ip2)  { 
ifp2=ifpl  <  1; 

theta=isi^*6.28318530717959/^fp2/ipl); 

wtemp=sin(0.5«theta); 

wpr  =  -2.0*wtemp*wtemp; 

wpi=8in(theta); 

wr=1.0; 

wi=0.0; 

for  (i3=l;i3<ifpl;i3+=ipl)  { 

for  (il=i3;il<i3+ipl-2aH-=2)  { 
for  (i2=il'42<ip3p2+=ifp2)  { 
kl=:i2; 
k2=kl+ifpl; 

tempr=wr*data[k2]-wi*data(k2+l]; 
tempi=wr*data[k2+l]+wi*data{k2]; 
datsi{k2] =data[k  1]  -  tempr; 
data  k2+l]=data(kl+l]-tempi; 
datakl]  +=  tempr; 
data[kl+l]  +=  tempi; 
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} 

} 

wr=(wtemp=wr)+wpr-wi*wpi+wr; 

wi=wi*wpr+wtemp*wpi+wi; 

} 

ifpl=ifp2; 

} 

nprev  ♦=  n; 

} 

} 


MAIN 

main(argc,argv) 
int  argc; 
char  ♦argvQ; 

{  FILE  *fin,  ♦fout; 

float  *output,*input,«truncjout; 
float  norm; 
float  *vector(); 
void  ♦free_vector(); 
char  name[30]; 

int  ij,  nn[i],  ndim,  isign,  new.order,  order,  imagejsize; 
if(argc  ^  5){ 

printf("!!!  The  command  line  should  be  !!!:\n\n  fft.trunc  infile  outfile  imsize 
order  \n\n''); 
exit(O); 

} 

image-size=atoi(  argv[3] ); 
order=atoi(argv{4]); 

sprintf(outname,  "X8Xd","netff t  .ng", order); 

up  dynamic  allocation****************^ 
input  =  vector(0,2*imagejize*imagejize-l); 
output  =  vector(0,imagejBize«image^ize— 1); 


/^,^t****************  Set  Up  Files  ************************/ 

if  ((fin=fopen(argv[l],"r"))  ==  NULL){ 
printf("I  can't  open  the  input  file"); 
exit(— 1); 

} 
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if  ((fout=fop€n(argv[2],"w"))  ==  NULL){ 
priiitf("I  can't  open  the  output  file"); 
exit(-l); 

} 

/*******^^***^^^L*Rea.d  File  ****i^**************i^ 

loopi(24>iinage-size*image_size— 1)  /*  initialize  array  to  zero  */ 

input[i]  =  0.0; 

loopi(imagejize*image^ize— 1)  /*re&d  data,  in  the  fourn  format  */ 
{  /♦  see  numerical  recipes  in  c  ♦/ 

f8canf(iin,  "HfXn",  &input[i*2]); 

} 


fclo6e(fin);  /*cIose  input  file  ♦/ 


/»  Initialization  parameters  for  FFT  */ 

nn[0]=imagejize;  /♦  size  of  mput  lAW  fournQ  ♦/ 
nn[l]=imagejsize; 

ndim=2;  /♦  two  dim  FFT  *f 

isign-l;  I*  FFT*! 

foum(input— 1,  nn-1,  2, 1); 


Find  Fourier  Magnitude 

j=0; 

for(i=0;i<(2>»imagejsize«iinage^ize-l);  i+=2) 

{ 

oatpat[i]=sqrt(  (double)  SQ(input[i])  + 

(double)  SQ(input[i+l])); 
j++; 

} 

norm=outpnt[0];  /*  d.c  component  used  for  normalization  ***/ 

printf("X4 .  Of  \n", norm); 

/*****  doUp******************************^ 

doflip(output4magej3ize);  /*  converts  fourn  format  to  human  format  *1 
printf(  "  X4 . 4f  \n  "  joutput  [8 1 28] ) ; 
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loopi(image^ize«im2ige  jsize)  { 
output[i]=  1000*output[i]/iiorm; 
fprintf(fout,  "Xf\n",  output[i]); 

} 

!*****  normalize  and  write  output  of  FFT  in  argv[2]  file  *****/ 

/*******  truncate  ******it******************^*^**t>^t*************************:^ 
/*  truncate  takes  fft(output)  of  size(image^ize)  and  truncates  the  FFT  to  **/ 

/*  order  specified  plus  d.c.  the  array  is  returned  in  trunc-out,  the  argv[2]*/ 

/*  is  used  as  a  header  when  truncate  writes  the  output  in  netfft.dat  */ 


if(order  0){ 

newjorder  =  2*order+l; 

truncjout  =  vector(04inagej5ize*image_size— 1); 

truncate(oatput  4mage^ize, order ,trunc-out,  argv[2] ); 


f*free.vector(  truncjout, 0,image^ize*imagejsize-l  );*/ 


fre€_vector(input  ,0,2«image  jsize*image^ize- 1 ); 
free.vector(output,04inagej5ize«imagejBize— 1); 
fclo8e(fout); 


} 


jacobi.c 
NAME:  jacobi.c 

DESCRIPTION:  Numerical  Recipies  routine  that  finds  eigenvectors  and  eigenvalues  of  a 
matrix 

SUBROUTINES  CALLED: 

WRITTEN  BY:  Numerical  Recipies  in  C 
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#include  <math.h> 


#define  ROTATE(a4j,k4)  g=a[i][j];h=a[lc][l);a[i][j]=g-s*(h+g*tau);\ 
a[k][l]=h+s*(g-h*tau); 

void  jacobi(a,n,d,v,nrot) 
float  **a,dD,**v; 
int  n,*nrot; 

{ 

int  j4q4p4; 

float  tresh, theta, tau,t,sm,84i,g,c,*b,^z,*vector(); 
void  nrerror(),free_vector(); 

b=vector(l,n); 
z=vector(l,n); 
for  (ip=l;ip<n4p++)  { 

for  (iq=l;iq<n4q++)  v[ip][iq]=0.0; 
v(ip]Iip]=1.0; 

} 

for  (ip=l;ip<n;ip++)  { 
b[ip]=d(ip]=a[ip][ip]; 
z[ip]=0.0; 

} 

♦nrot=0; 

for  (i=:l;i<50;i++)  { 

8m=0.0; 

for  (ip=l;ip<n-14p++)  { 
for  (iq=ip+l;iq<n;iq++) 

6in  +=  fab8(a[ip][iq])’, 

if  (8m  ==  0.0)  { 
free-vector(z ,  I  ,n) ; 
free_vector(b,l  ,n); 
return; 

} 

if  (i  <  4) 

tresh=0.2*8m/(n*n); 

el8e 

tre8h=0.0; 

for  (ip=l;ip<n-l;ip++)  { 
for  (iq=ip+l;iq<n;iq++)  { 
g=  100.0*fab8(a[ipl[iq]); 
if  (i  >  4  &&  fab8(d[ipj)+g  ==  fabs(d[ip]) 

&&  fab8(d[iq[])+g  ==  fab8(d[iq])) 
a[ip][iq]=0.0; 

el8€  if  (fab8(a[ip][iq])  >  treeh)  { 
h=d[iq]-d[ipj; 
if  (fab8(h)+g  ==  fab8(h)) 
t=(a[ip]liq])/h; 
el8e  { 

theta=0.5*h/(a(ip][iq]); 
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t=1.0/(fabs(theta)+6qrt(1.0+theta*theta)); 
if  (theta  <  0.0)  t  =  — t; 

} 

c=1.0/8qrt(l+t*t); 

8=t*c; 

tau=8/(1.0+c); 
h=t*a[ip][iq]; 
z  ip]  -=  h; 
z  iqj  +=  h; 
d[ip  -=  h; 
d  iq  +=  h; 
a  ip]liqj=0.0; 
for  0=lii^‘P~llj++)  { 

ROTATE(  aj  4p  J  m) 

} 

for  (j=ip+llj<iq-lu++)  { 

ROTATE(a4pJj4q) 

} 

for  U=iq+llj^nU++)  { 

iiOTATE(a4pJ4qj) 

} 

for  (j=lu<n-J++)  { 

ROTATE(v44pJ4q) 

} 

++(*nrot); 

} 

} 

} 

for  (ip=l4p^n;ip-H+)  { 
b(ip]  +=  z[ip]; 
d[ipl=b[ip]; 
z[ip]=0.0; 

} 

} 

nrerror("Too  many  iterations  in  routine  JACOBI"); 

} 

#undef  ROTATE 


eigsrt.c 

^*#«*«**************4t*«*«««««****««*«***««*4t«***«*«*«4r«**4>«***4t*#**«#4<#4>4'*t 

NAME:  eigsrt.c 

DESCRIPTION:  Numerical  Recipies  routine  that  sorts  eigenvalues  and  vectors  into  decr- 
reasing  order. 

SUBROUTINES  CALLED: 

WRITTEN  BY:  Numerical  Recipies  in  C 


***************************^1****************************************’^ 
void  «gsrt(d,v,n) 


float  d[],**v; 
int  n; 

{ 

int  kJ4; 
float  p; 


for  (i=l;i<n;i++)  { 
p=d[k=i]; 

for  (j=i+lu<n'J++) 
if  (dli]  >  p)  p=d[k=j]; 
if(k#i){ 
d[k]=d[i]; 
d[ii=p; 

for  (j=lu<nu++)  { 


P=vu  I 

vDl[il=v 

vS][k]=p; 


fourn.c 

/♦♦*♦♦*♦***♦♦♦*♦♦♦*♦***♦♦♦*♦♦*♦*♦♦♦♦♦♦*♦♦**♦***♦**♦♦*♦♦*♦*****♦*♦****♦♦♦♦** 

NAME:  foura.c 

DESCRIPTION;  Numerical  Recipies  multi  dimeasional  FFT  routine. 

Requires  a  complex  column  vector  as  follows: 

/  real  a(l)/ 

/  complex  a(l )/ 

/  real  a(2)/ 

/  complex  a(2)/ 

/etc/ 

SUBROUTINES  CALLED: 

WRITTEN  BY:  Numerical  Recipies  in  C 

#mclude  <math.h> 

#define  SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr 

void  fourn(data,nn,ndim4sign) 

float  dataQ; 

int  nnQ,ndim4sign; 

{ 

int  il424342rev43rev4pl4p24p34fpf4fp2; 
intibit4dim,kl4c2,n,nprev  ,nrem  ,ntot ; 
float  tempi, tempr; 
double  theta, wi,wpi,wpr,wr,wtemp; 
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ntot=l; 

for  (i<lim=l;idiin<ndim;idim++) 
ntot  *=  nn[idim]; 
nprev=l; 

for  (idiin=ndiin;idim>l;i<lim — )  { 
n=nn[idim]; 
iirem=ntot/(n*nprev); 
ipl=nprev  <  1; 
ip2=ipl«n; 
ip3=ip2'»nrem; 
i2rev=l; 

for  (i2=l;i2<ip2;i2+=ipl)  { 
if  (i2  <  i2rev)  { 

for  (il=i2;il<i2+ipl-2;il+=2)  { 
for  (i3=il;i3<ip3;i3+=ip2)  { 
i3rev=i2rev+i3— i2; 

SWAP(data  i3],data[i3rev]); 
SWAP(data  i3+l],data[i3rev+l]); 


ibit=ip2  >  1; 

while  (ibit  >  ipl  &&  i2rev  >  ibit)  { 
i2rev  -=  ibit; 
ibit  >=  1; 

} 

i2rev  +=  ibit; 

} 

ifpl=ipl; 

while  (ifpl  <  ip2)  { 
ifp2=ifpl  <  1; 

theta=i8ip*6.28318530717959/Cifp2/^l); 

wteinp=8in(0.5«theta); 

wpr  =  — 2.0*wtemp*wtemp; 

wpi=8in(theta); 

wr=1.0; 

wi=0.0; 

for  (i3=l;i3<ifpl;i3+=ipli  { 

for  (il=i3;il<i3+ipl-2;il+=2)  { 
for  (i2=il;i2<ip3;i2H-=ifp2)  { 
kl=i2; 
k2=kl+ifpl; 

tempr = wr '*data[k2]  -  wi+data[k2+ 1]; 
tempi=wr«data[k2+l]+wi*data[k2]; 
data  k2]=data[kl]— tempr; 
data  k2+l]=data[kl'|-l]— tempi; 
data  kl]  +=  tempr; 
data[kl+l]  +=  tempi; 

} 

wr=(wtemp=wr)*wpr— wi*wpi+wr; 
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wi=wi*wpr+wtemp*wpi+wi; 

} 

ifpl=ifp2; 

} 

nprev  *=  n; 

} 

} 

#uiidef  SWAP 


nrutil.h 

float  ♦vector(); 
float  ’•■*matrix(); 
float  **convert  jnatrix(); 
double  *dvector(); 
double  **dinatrix(); 
int  ♦ivector(); 
int  **imatrix(); 
float  ♦^submatrixO; 
void  free_vector(); 
void  free_dvector(); 
void  free  jvector(); 
void  free-matrix(); 
void  free-dmatrix(); 
void  freeJmatrbc(); 
void  free-8ubmatrix(); 
void  free_convert_matrix(); 
void  nrerror(); 
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